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Chapter 8. FALLACY OF APPLYING STAGGERED GRIDS 
IN NUMERICAL CALCULATIONS 

§1. Irreducible errors in schemes built upon staggered grids 

We'll show irreducible errors on Harlow-Welch scheme. Harlow- 
Welch scheme III uses in two-dimensional problems three staggered 
in relation to each other net domains, in three-dimensional problems 
- four. Let's limit ourselves to consideration of two-dimensional 
problem for Navier equations of incompressible fluid 

du du 2 duv 1 dp ,d 2 u SV 
ot ox dy p ox ox oy 



dv dvu dv I dp ,d v d v 

— + + + — = v( — j + — 

dt dx dy p dy dx dy 



+ — + — + -^f = K^T + ^r)' (2) 



du dv 

— + — = (3) 

dx dy 

with initial u\ t=Q = d u (x, y),v\ t=Q = d v (x, y), 

and boundary conditions U\ s = <p u (x,y,t),V\ s = <p v (x,y,t) . 

In rectangular domain Q = [0<x<X, 0<y<Y] with 

boundary S nodes of the major grid 

n h ={x i =ih x ,i = 0,N x ;y j =jh y , j = 0,N y } 
are such that 



{x =0,x Nx =X,j = 0,N y },{y =0,y Ny =Y,i = 0,N x } 

are boundary nodes lying on physical boundary of the domain. Two 
additional grids Q M , Q v are staggered relative to this main grid: 

U u ={x l =(i + -)h x ,i = -l,N x ;y j = jh j = 0,N y }, 

>'+2 2 ■ 

1 



& v ={x i =ih x ,i = 0,N x ;y 1 =(j + -)h y , j = -l,N y }, 

J+ 2 L 
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Sets of internal nodes of these three grids are denoted accordingly 

& h = k = *K l = ! ' N * - 1 ; y.i = ; -/ v, J = l > N y -v> 

a=U. l =d + bh x ,i = 0,N x -l;y j =jh y , j = l,N y -l}, 

' + 2 2 

1 



Q v ={ Xi =ih x ,i = l,N x -l;y l =(j + -)h y , j = 0,N y -l}, 

J+ 2 2 ' 

Timegrid Q r ={t n =nr, n = 0,l,--,N T } . 

Used are index notations of grid functions 

fZ=f<,x t ,y p t n ),f\ =f(x x ,y r t n ), 

i±-j i±- 

2 2 

f\=f(x i ,y x ,t n ),f\ ,=/(* !,> lO. 

2 2 2 2 2 2 

The Harlow-Welch scheme uses three grids, its own grid domain is 
built for each equation. Therefore equation (1) approximated in 

nodes of grid Q u : 

n+\ n n n n n 

U ! — U , w2 „2 M 1 1 V 1 1 _M 1 1 V 1 1 

H — / H — / ll- , • — li- H — H — H — /H — H — / H — / 

2 2 , '+lj !/ , 2 2 2 2 2 2 2 2 , 



n „ U 3 -Ml j + W j Mj -2« , +« j 

, p '+^~^ =v( 'V IV |V | 'V' +1 V 'V'"\ 

/ = 0,...,JV X -1, j=l,...,N y -l, (4) 

equation (2) is approximated in nodes of grid Q v : 

n+\ n n n n n 

V .. l-V.. l «. 1. lV !.!-«!. iV ! ! „ 2 n2 

»+- </+- (+-;+- i+-/+- i — ;+- i — ;+- V-. , — V- 

2 2 , 2 2 2 2 2 2 2 2 , J/+ 1 V , 



*, *, 
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Fg+1 Fg _ , J 2 ' 2 J 2 , J 2 y 2 J 2 \ 

i = l,...,N x -l, j = 0,...,N y -l, (5) 

continuity equation (3) is approximated in nodes of the main grid 

a 



'h ■ 



u n+ \ -u n+ l v" + \ -v n+ \ 

i+ 2 j '-2 j , ii+ 2 iJ -2 



+ ^ *- = 0,i = 0,N x J = 0,N y (6) 



*, hy 



With values i = 0, j = l,...,N — 1, (6) includes N —1 values 
u n f and when i =N x ,j = l,...,N -1 (6) includes N y -1 of 

~2 J 

,.«+! — 

values M l . of unknown functions beyond physical domain Q . 

' + 2 } 

Similarly, when values j = 0,i =l,...,N x —1 , (6) includes A^ —1 
of values v l and when j = N , i = 1,...,N X — 1 , (6) includes 

i — N,. 

2 J 

«+l 

N x —1 of values v . l of unknown functions beyond physical 

y 2 

domain Q . Along with this unknown u j are included into (4) 

-2 J 
,.1+1 

when i = 0, / = 1.....A 7 ,, — 1, and u l are included into (4) when 

2 

i=N x -l, j = l,...,N y -l. 

Here arises a natural question: how to define them? Equations (6) are 
not suitable for this purpose. According to the general concept, these 
equations are necessary for computation of pressure 
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P ;,i = 0,...,N x , j=0,...,N y 

(nodes on rectangle vertices are not included in here). 
From equations (4) and (5) singled out are 



„+l Pi+lj Pij 



V fix V 



+ R% ,i = 0,...,N x -l,j=l,...,N y -l, (7) 



v n+ \ =-r Fii+ \ Fij +Q n , ,i = \,...,N x -l,j = 0,...,N y -1 (8) 

iJ + 2 Phy iJ+ 2 

and are substituted into continuity equation (6) in nodes with 
numbers i = 1, . . , ,N X — 1, 7 = 1, . . . , N — 1 . In the result obtained is 
a system of difference equations for pressure 

p;«j -w + pt-u pu-w+pu 



+■ 



ph 2 x ph 2 y 

R\ -R\ Q\-Q n , 

'■1 — J ' — i b'h — ') — 

2 + 2 ^ 2 - , (9) 



i=l,...,N x -l, j=l,...,N y -l 

Boundary conditions for pressure are derived from the same 
equations (6), but considered in boundary nodes 

i = 0,j = l,...,N y -l;i = 0,j = l,...,N y -l; 

j = 0,i = l,...,N x -l;j=N y ,i = l,...,N x -l 

Let's limit to boundary nodes i = 0, j = l,...,N — 1, to show 

arisen problem of defining u { . 

Substituting when i = 0, j = l,...,N — 1 into continuity equation (6) 
values (7), we obtain differential boundary conditions for pressure 
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n n 



(-r Plj Poj +r; )-u"f v" + \-v" + \ 

rjfi -j — j oy+- oy — 

^— ~ ~+— *; 2 - = 0,j = l,N -1, (10) 

K h y 

where due to boundary conditions for velocity components it is 
provided 

v " + ji = <Pv(*o>)\ + nO, j = l,...,N y -l 

As can be seen, boundary condition (10) includes unknown 
U j , j = l,...,N — 1 . Similar will be entry of unknown quantities 

~2 J 

into boundary conditions for pressure 

" ^ j = l,...,N y -l;v+l,v ^,i = l,...,N x -U 

' 2 2 ' 2 



§2. On additional boundary conditions and on the problem of 
one-valued solvability 

So, for defining unknown quantities 

u n f,u n+ \ ,j = l,...,N y -l,v n+ l,v n+i „i = \,...,N x -\ 

— ; N r +-j i — iN x +- 

2 * 2 2 " 2 

required are additional boundary conditions for calculations, which 
contradicts the initial setting of initially-boundary problem for 
equations (1),(2),(3), because there are enough available boundary 

conditions u\ s = <p u (x, y, t),v\ s =<p v (x, y, t) . 

For example, Poisson equation with boundary condition of Dirichlet 

a 2 ® s 2 ® , 

ox oy 
has a unique solution, while in applying staggered grids for defining 

®";\®" +1 , J=l,...,N y -l, 0» + j,0" +1 lt i = l...,N x -l, 

— j N x +-j i — iN v +- 

2 2 2 2 

Required is additionally boundary condition of A®L = (// type, 
where, furthermore, arises a problem of defining type of differential 
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operator A and function iff = l//(x, y) on boundary of domain S . 
And right away arises a question: is there such unique solution of 
Poisson equation with two boundary conditions <X>L = (p and 

AOL = (//, which would match with solution of Poisson equation 

with one boundary condition OL= (p! 

Literature 
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Chapter 9. TECHNOLOGIES OF EFFECTIVE SCHEMES 
BUILDING ON A SINGLE GRID 

Families of schemes, using only one main grid Q. h for equations 

Navier-Stokes were developed and theoretically justified in 
monograph of the author III. 

Let's stop only on some of them on example of a two-dimensional 
problem. Calculation of disk longitudinal flow, results of which 
confirmed fallacy of Prandtl equations, was carried out with the use 
of one of these schemes. 

§1. Formulation of the problem longitudinal flow past a plate 

For calculation of two-dimensional flow past a plate, let's use 
closed system of equations from Chapter 3: 
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pV— + > v,.^ + > — ^- > '— ] + —-- 

= pF i + J uAv i , i= 1,2,3, 

p — -+y (v. — L + v', — L -)H — > — ] + — 

■* y 7 ■* y ( 



(i) 

(2) 



= pF' i+J uAv\, 



1,2, 



Yax, 



,= ?>',- 



,o=V? 



(3) 

(4) 

(5) 



Via =( Pi^i\t=0= V i^i 
for simulating of both laminar ( v\ = 0,i = 1,2,3 ) as well as turbulent 
( v\ ^ 0, i = 1,2,3 ) modes (interleaved mode is accounted for 
automatically) of flow past a plate and other objects. In notations 
v l =u,v 2 =v,p = p,F l =F x ,F 2 =F y ,v\ = u',v' 2 = v' system will 
take the form of 



,du du du du a du'v' t° d ,du a du'v' 



p{ — + u — + v — + + 



- + - 



dt dx dy dx dy 2 dt dx dy 



)} + 



dx 



P F *+ju( 



d 2 u d 2 u 

dx 2 dy 2 



), 



(6) 



,dv dv dv du'v 1 dv t d ,duv dv . . dp 

p{— + u — + v — + + ( + )} + — 

dt dx dy dx dy 2 dt dx dy dy 

( d 2 v dV 



(7) 
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^ + ^=0, (8) 

dx dy 

,du [ du' du [ .du .du t° d ,du a du'v'., dp 1 
p{ — + u — + v — + u — + v — + ( + )} + — = 

dt dx dy dx dy 2 dt dx dy dx 



av d 2 u\ 

— T + — ; 
dx dy' 



= P f x +M(^^ + ^tt)' (9) 



r dv' dv' dv [ ,dv ,dv f d ,du'V 5v' 2 xnl dp' 
p{ — + u — + v — + u — + v — + ( + )]}+— = 

dt dx dy dx dy 2 dt dx dy dy 

, ( d 2 V d 2 v\ 
= pF + /j(— T + — -), (10) 

ox dy 

^ + ^ = 0, (11) 

dx dy 

Equations (6),(7),(8) transform into Navier incompressible fluids 
equations when u'=0,v'=0,V(x,y,t) . 
At the zero time t = fluid is in state of rest in laminar flow: 
u(x, y,0) = 0, v(x, y,0) = 0,0<x<L + 21,0 <y<H, 

at some moment of time t = t * by flow obtained for this moment of 

n 

time by numerical method, imposed is a "random" disturbance, i.e. in 
the grid nodes at this moment are assumed 

* * 

u?j =u%+g p vl =v$+{ ij ,i=l,..,N x -l,j=l,...,N y -l, 

where g ij =g(x i ,y j ) and <^ i j=^(x i ,y j ) close to zero random 

functions, simulating initial perturbation in a flow (they can be set 
with the help of random numbers generator). 

Boundary conditions: incident flow has the following velocities 
x = 0,w(0, y,t) = U x exp(-Z?70,v(0, y,t) = 0, 

u'(0,y,t) = c]{y,t) -random function, v'(0,y,t)=0, (12) 

0<y <H,b'= const > 0, \^(y,t)\«U a0 eet., 
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on the plate - condition of adhesion and impermeability 

y = Q,u(x,0,t)=0,u'(x,0,t) = 

v(x,0,t)=0,v'(x,0,t) = 0,l<x<L + l, (13) 

at the flow outlet "soft boundary conditions" are set 

a 2 a 2 

BML + 2l,y,,) BW(L + 2l,y,,) 

OX ox 

in the upper part of flow, extremum conditions 

H du{x,H,t) _ du'(x,H,t) _ Q 

dy dy 

^^H,t) =0 dv'(x,H,t) =0Q ^ x ^ L + 2l (15) 

dy dy 

Before beginning of the plate and after the plate symmetric flow- 
around conditions are set 

du(x,0,t) du\x,0,t) 
y = 0, = 0,v(*,0,0 = 0, = 0,v 0,0,0 = 0, 

dy dy 

0<x<l,L + l<x<L + 2l 

Transition to dimensionless variables 

For transition to dimensionless variables selected as scale values 

are: for velocities U x , for linear dimensions - plate length L , for 

time L/U m , for pressure pU x . In dimensionless quantities, 

£=tU x /L,£=x/L, 

fZ=y/L,£=u/U x ,€=v/U x ,£ = p/(pU 2 J,Dg=t'VJL, 

F x = F x / g.F y = F y I g.Fr = U 2 l{gL) — Froude number, 

Re = pU x L/ jU system (6), (7), (8) (signs are omitted) has the form: 

du du du du' 2 du'v' Dg d ,du' 2 du'v\ dp 

— + u — + V — + + -— ( + ) + — = 

dt dx dy dx dy 2 dt dx dy dx 



175 



F 1 / d 2 u 5 V 

= — - + — ( — y + — r)» ( 16 ) 

Fr Re dx dy 

dv dv dv du'V dv' 2 Dg d ,du'v' dv' 2 ^ dp 

— + u — + v — + + — ( + ) + — = 

dt dx dy dx dy 2 dt dx dy dy 

F y 1 d 2 v 5 V 

= — + ( — v + — 5"). (17) 

Fr Re dx dy 

8 l + ^=0 (18) 

dx dy 

Equations for disturbances (pulsations) (9),(10),(1 1) are also 
transformed into similar dimensionless form: 

du' du' du' ,du ,du Dg d y du' 2 du'v\ dp' 

— + u — + v — + u — + v — + ^-— ( + ) + — = 

dt dx dy dx dy 2 dt dx dy dx 

F' 1 ,5V d 2 u\ 
= — - + — ( — t + — r)» ( 19 ) 

Fr Re dx dy 

dv' dv' dv' .dv .dv Dg d ,du'V dv' 2 . dp' 

— + u — + v — + u — + v — + ^- — ( + ) + — = 

dt dx dy dx dy 2 dt dx dy dy 

Fy 1 ,aV d 2 v\ 
= 1 ( — t + — r)» ( 2 °) 

Fr Re dx dy 

* + ^ = 0, (21, 

Boundary conditions (12)-(15) 

x = 0, w(0, y,t) = exp(-b°/t), 

v(0,y,t)=0,0<y<H/L,b° =const >0; 

dx 2 
S\a + 2,/L,y, l) 

OX 
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}' = 0, u(x,0,t)=0, v(x,0,t) = 0,UL<x<l + l/L; 

du(x,H/L,t) 

y = H I L, = 0, 

dy 

g^^M = 0,0<x<l + 2//L; 

dy 

du{x,Q,t) 

dy 

v(x,0,t)=0,0<x<l/L,\ + l/L<x<\ + 2l/L 
Similar dimensionless boundary conditions occur for pulsation 

equations as well. 
Initially-boundary problem for (16)-(21) is solved theoretically by 

justified difference methods Jakupov K.B. Ill), detailed technique of 

which is given below. 

§2. Technologies of building difference schemes in vicious fluid 
dynamic equations on a single grid. Semi-implicit difference 

scheme 

For stating principles of building and algorithm of realizing 
difference schemes for vicious fluid dynamics equations on a single 
grid, the most suitable is the family of parametric semi-implicit 

schemes from /l/( a = 0, /? =1), which is represented by example of 

equations (16), (17), (18). In dimensionless domain 0<x<l+2ll L, 
0<y<H/L grid U h ={x i J = 0,...,N x ;y j J = 0,...,N y } is 
introduced, and time grid Q^o = f^,n = 0,l,..., N T with pitches 

K = x i ~ x i-i >°> i = l->N x ,h yj = y } - y H >0,j = l,...,N y , 

x =0,y =0,x Nx =l+2l/L,y Ny =H/L,T n =t n -t n _ l >0,n = l,...,N z 

Grid pitches near streamlined surface must be of the following order 

h « 1 / vRe . Monotonic scheme Chapter 6 approximating 
convection members with 2 nd order of accuracy is applied (see 
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Jakupov K.B.I2I). Scheme is semi-implicit in time because pressure 
and continuity equation are taken on the upper time layer: 

n+l n \ n \ . n n \ n \ \ n i . n n \ n \ 

B B , H i n , i i n , ij i] n , i i] n , 

-H — L u-+— —u+ — -u-+— —u + 



2 2 2 2 



v 



+ {(W )j+(KV)j — [(« )j+(MV)y-(M )j ~{UV)~ ]/T n } + 



Pi+lj Pij _ F x 2/4„ , „ „. , 2 /Vy« 



n+l _ n+l 

+ M ; * =if-+T^-(«:-o+ "^ «-<)],(22) 
^,+i *> ^+i+^« Vi+\- 



n+l n \ n \ , n n \ n \ \ n i , n n \ n \ 

V.. -V ;; W,.. +W-. W,, - «,.; V,, +V,, V.. - V.. 

9 9 i 9 9 1 i 9 9 1 i 9 9 » i 9 9 n , 

— -H -V- H — — V H -V-+ — — V + 



^n+\ 



,i2 \n 6 r/ i,.i\h I /,.'2\n / t„.i\»— 1 / i2\rc— In 



+ {( M 'v')'~+(v' z );-^[(M'v');+(v' 2 );-( M 'v')r-(v' 2 )r 1 ]/r„} + 

+ — ' — = — + — (v -v =) + — (v -v-)J, (23) 

h Fr h +h h +h 

i = l,...,N x -l,j = l,...,N y -l; 

u n+{ -u n+ }. v H+1 -v" + ! 

V >-l y , !/ 9-1 



+ - , ' =0, i = U.,N x ,j = l,...,N y ; n = 0,N T , (24) 

where coefficients in dissipative members are equal 

\ u ;\+u" \u"\-u" lv"l+v" lv-1-v" 

91 = 1 + Re(^ n xi + +^ l! -h KM + -? ?- A . + -? ?- A . +1 ), 

4 w 4 -" +1 4 w 4 w+1 

1 I v" I +v" I v" I -v" 

Re 4 4 n yj+i' /JX ' 

1 Im" I+m" Im" I-m" 

r"vyu V -p, . *xt A Xl+l' 

Re 4 4 
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On uniform mesh h yj+l =h yj =h y ,\/j,h xM =h xi =h x ,\/i these 
coefficients are reduced 

1 \v-Ah \u"\h r \v"\h v 

Re 2 2 2 



1 \u"\h \u"\h r \v"\h 

ju" =[— + -^] /[l + Re(^^ + ~^— 

vyu Re 2 2 2 

With values 

Re ' Mvyu Re 
the scheme has the 1 st order of accuracy. 
Initial conditions are set on grid 

a h =$=l,...,N x -l;y p j = l,...,N y -l: 



r^uxv -pi _ ' r^vyu 



ul=0,v°=0,i = \,N x -\,j = \,N y -\ (25) 

Boundary conditions (18)-(21) are set in boundary nodes (the below 
formulae are given on uniform grid): 

S h ={i = 0,i = N x J = W y ;j = 0,j = N y ,i=W x },O h =n h [jS h ; 

b" 
x = 0, u n Q f = e '"♦' , VoJ 1 =0,j = 0,N y -l, (26) 

(27) 



1 , r\i I j n+1 r\ n+1 n+1 

X — 1 t Li I Lt, U N j — LU N _ 1 • — U N _2j •> 



(28) 



n+i r\ n+i n+i c\ ^ • • at 

v nj = 2v N t -ij- v N x -2j$ ^J^N y ; 

on plate w" +1 = 0, v" +1 = 0, mk < i < mkk, 

before the beginning and after the end of the plate njiacTHHbi 

y = 0, u" +1 = (Auf - uf ) 1 3, v,; +1 = 0,1 < i < mk - 1, mkk + 1 < i < N x , 

TJ 

n+1 / a n+1 n+1 \ / n 

y= — - U iN y = (Hv v -1 - U iN-l) I 3 > 

L (29) 

v:;;=(4v:;;_ i - v ^_ 2 )/3j=o,N x , 

y y y 

Scheme for pulsation equations has the following form: 
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tM+1 m i n i . n n \ n \ \ n \ . n n \ n \ 

U ,.; ~U ,, M,.. +«;; W, - «;; V;; +V;; V ;; ~ V;; 

-* ' L + ^ ' L U'1 + ^ '-14'* + -* ' L u'" + ^ ?- M '» + 

z n+l 2 2 2 J 2 ' 

+ u ij u~ x +v ij u-+—[(u )i+(«v) r (« )j -(«v) 5 ]/t„ + 

m+1 — n'" +1 E"' 7 l ' 2/i" 

+ ^. + U gjj = ^ + fAgv ( M '"- M '")+ ^>" «- M '")], 

v'" +1 -v'" lw"l+w" m"-Im"I I v" I +v" v"-lv"l 

g g , g V . .in . g g in i/ i/ in , !/ 'J in , 

r , 2 * 2 J 2 y " 2 



i «+i 



.» , D s 



+ M # V S +V # Vy+— [(»V)j+(V )j-(KV)i ~(V )y ]/T„ + 

+ r v+ i Pa = ^ + fAov ( V ';- V ^)+ "va ( v ' a -v'" ? )], 
i = l,... ,^-l,7=l,...,AT y -1; 

g - '~ lJ + iJ - ij ~ l =0, i=l,...,N x J = l,...,N- n = 0,N 

§3. Semi-implicit scheme realization algorithm 

Scheme realization algorithm (22)-(29) considerably simplifies, 
if notations of algebraic sums of convection and dissipative members 
are introduced into (22), (23): 



u " Km+K Km K, h yj + i+K K + i K 
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i«;i«;«;-«^ «;-i«;i«^-«; 



2 ^ 2 ^ + i 

|y.. | +Vj . «..-«.,. Vg-lv^lM^,-^. 

2 ^ 2 ^. +1 

2// v v",.-v" v"-v",. 2//" v",-v" v"-v", 
" Km+K Km K Km+Kj Km K 

\u;\ + u; v ;-v: Aj *;-i«;iv^-v; 



(31) 



2 K 2 /z n+1 

| v "| +v " v "_ v " v "-lv"lv" ,-v" 

2 ^ 2 ^. +] 

Schemes (22), (23) are reduced to a standard form 

n+l n „ B+1 »,"+! 

v " +1 - v n D n+l — D n+l 

^^L = Q:„- P -\ P " ■ (33) 

T n+\ yj+i 



where 



-^[( M ,2 );-+( M v);-( M ,2 )r 1 -( M , v)r 1 ]/r„}, (34) 

- *y [(«' v ) ~ + ( v 2 ); - (u' V ' y.- 1 - ( v' 2 );-' ] / r„ } (35) 

Both parts (32), (33) are multiplied by r„ +1 and are represented in 
the form 
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n+1 n " +1 

n+1 n , _ 7T« ri+lj rij ,.„ 

W ? y =W/,-+^+lG^-Vi 7 — • ( 36 ) 

^i-i+i 

n+1 n+1 

n+1 n . _ 7T« rii+\ Pij ,„-, 

V y = V ij + ^n+lGvy ~ ^„+l ^T ~ ( 37 ) 

In (36) and (37) let's introduce with the purpose of shortness, 

Kj = K +r n+l Q u n ij , F^=V1 +r^O; - (38) 
and represent (36),(37) in even shorter form 

C" +1 - D n+l 



r,+l; J 


"V 


^,+1 

n+1 


Pij 



n+1 77. n -r J/+1 r l] 

v a = F vij - r n+ i , (40) 

Absolutely similar (30) - (40) expressions are obtained for the 
scheme of pulsation equations as well (19), (20), (21): 

u m+1 -u m - P m ll-P m+1 

_g _g_ _ ft <n r ;+ly f ij 

t j h 

4 n+1 n xi+\ 

v m+1 -v m — p m+ \-p m+l 

<j JL — ft m y+1 y 



where 



t ~ vlJ h 

'n+1 "yy+l 



— F 

— O ~-\ — -—uu~—vu~ — 

mj ^ mj p U x y ij y 

- — [(w )~+(uv)~-(u )- -{uv)~ ]/r n , 

Q<» = Q m .. + ^-u m V" -V'" V" - 

& r/ 1 i \n 1/ i2 \ n / t i\ n— 1 / i2 \n— l-i / _ 
y[("V);+(v )--( MV )- -(V )j ]/T a 



Having denoted 
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F uij = U ij +T n+ lQ uij . F vij = v a +T n+ lQ vij . 
we'll make a standard representation 

m+l _„m+l 
m+l _ 77m _ £_ i+l 7 ^ g 

M y — r uij Z n+\ h 

m+l m+l 

in+1 _ 77 1« _ i' y+1 i' y 



'7 vij ""n+1 , 

There is a fundamental principle of realization of law of mass 
conservation via pressure function. By this principle pressure field 

p", + is determined from requirement, that for those standing in (39) 

and (40) w" + , v" + observed would be continuity equation (24). With 

this purpose (39), (40) are substituted into (24), and in the result 
obtained is the closed system of difference equations for 

pressure p"* . Substitution is performed gradually. 

§4. Technology of obtaining difference equations for pressure is 
explained for boundary conditions 

u\=(p u {x,y,t),v\ s =(p v {x,y,t) 

l nd scheme 

1". Substitution of (39) and (40) into (24) in nodes with numbers 

2 < i < N x -1,2 < / < N - 1 gives a system of linear equations 



1 n+1 n+1 n+1 n+1 

— 4f"- t Pi+li ~ Pij \-€F n -r Pii ~ Pi ~ lj ) + 

n xi "iM "xi 

n+1 «+l n+1 n+1 



1 r\ v\ 

"'" t ^ i">' T n+1 , / " r vi/-l ^n+1 , -* J" u ' 

2<i<N x -l,2<j<N y -l (41) 
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2°. Substitution of (39) and (40) into (24) in nodes with indices 
1=1; 2 < j < N — 1 . In these nodes equation (24) has the following 
form: 

n+l n+l n+1 n+1 

U ^^ + V -^^± = 0,2<j<N -1, (42) 

where boundary condition w . = <p u0j is set. After substitution the 
result is: 

n+1 n+1 -. 

— <F ulj - T n+l )-u 0j +— *F v ij ~ (43) 

n xl n x2 rl yj 

n+l «"+! tt +l n +l 

Pi;+i ~ Pi y s y ffl Py_~Piy-K 1 n o <r ; <r a/ 1 

- Vi : ) - i A-ii-i - Vi : ) yO,2<j<N - 1 

3". Similar substitution (39) and (40) into (24) is performed in 
nodes with indexes j = 1; 2 < i < N x — 1, where (24) has the form 

n+1 n+l n+1 n+1 

W , — U ,, V-, — V- n 

_d di + ^l 15_=0, 2<i<N -1, (44) 

n+1 _ n+1 

here, in accordance with boundary condition given is V i0 — (p vi0 . 
In the result 

1 n n+l - n" +1 n" +1 _ n" +1 

, X- 1 Mil 'n+1 , / V am '„+i J/^ 

^xi ^.ri+1 "xi 

+ 7 L {(^i-r„ + /'\~ A " +1 )-v i ; +1 }-0,2</<iV j -l (45) 

4". Substitution of (39) and (40) into (24) in corner node with 
indexes i = 1; j = 1; where (24) has the form 

n+1 n+1 n+l n+l 

M ll — M 01 . V ll ~ V 10 _Q 



^1 h yl 



results in expression: 
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1 r," +1 _ n" +1 

1 4f" -t Pn Pn 

X 1 nil 'n 



K * "" " +l K 2 )-< 3 

i tt+1 «+l 

+ - -l^i- r^^-p^-)^ 1 30 (46) 

In (41), (43), (45), (46) substitution in internal nodes Cl h was 
carried out, the resulted system of linear equations for p" + is open- 
ended. For closing this system, engaged is continuity equation (24) in 
boundary nodes y h with indices i = N x , l<j<N -1: j = N , 

2<i<N -1. 

5". In boundary nodes i = N x , \<j<N -1 continuity equation 
(24) has the form 

,,"+! _,." +1 n+1 n+1 

" J ^/^ rV =0, l<;<iV-l, (47) 

where by boundary conditions the following considered known 
n+1 „n+l n+1 „fl+l n+1 ,«n+l ■ a at i 

«*u = W%- =fl*y. v-i =^-i .; =o,tf, -i, 

n+1 

accordingly, only (39) is substituted instead of W^ _^ • . As the 
result obtained is 

- — l^j-^-ij-^+i ; )}+ ; = °> (48) 

l<j<iV v -l 
d". In boundary nodes j = N , \<i<N x -\ continuity equation 

(24) has the form 

n+1 n+1 n+1 n+1 

M ;W V ~ u i-\N v v iN v ~ v iN v -l 

—?- y - + ^ — = 0, l<i<N x -l, 

K h yN y 

where by boundary conditions the following considered known 
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n+l n+l «+l n+l n+l n+l • r\ »r 

U iN y = <P*N, » V «, = <PvW y . «,--W, = fto-W, , I = 0,. . .,Ar x , 

n+l 

therefore (40) is substituted instead of ^j/v -i . As the result 

y 

obtained is 

n+l n+l n+l n+l 

-^ L +— <-(C,-i-^n+i , )} = 0,1<»<JV,-1(49) 

Difference equations (48) and (49) are natural boundary conditions 
for pressure. System (41),(43),(45),(46),(48),(49) is the closed 
system, which is solved by iterative method, for applying of which it 
is convenient to represent equations (41),(43),(45),(46) in the form of 

one formula using signature function sign 4 '■ 

sign 4* =1 for x > , sign 4 =0 for x = , sign 4 j= —1 for x <0, 



J_ Jtpn _ Pi+lj-Py v 1 + Sign %,-l.5^ 

7 ^ uii n+l i ' ft * ui-lj 

K Km 2 

n+l _ n+l 

_ T Ll LtLL\ l-sign(i-l.5) B+1 -j 1 - 

hi 2 0j J A H . ^ " J 

_ r Py + i~Pi; x 1 + Slgn f-1.5^ „ _ Py ~ Pq-i 

-Vi ^ )- 2 (**-, ^i ^ )" 

_ l- J »ffl,(j-l.S) v . +1 3o,l</<iV x -l,l< 7 -<^-l (50) 



System (50) is solved in relation to p" + with boundary conditions 
(48), (49). 

2 nd scheme with forwardwise difference derivatives in continuity 
equation ( a =1,/? =0 /l/): 

u n+l — u n+l v " +1 -v" +1 

^+17 ^ + ^/±i '^ = 0,0<i<N x -l,0<j<N -1 (51) 
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In this case pressure gradients in equations (16) and (17), 111, 
should be approximated by backward differences, 

n+l n+l 

K 

n+l n+l 

vT l =K v -r J" ~ Pii ~ l (53) 

K 

In the result of similar 1°, 2°, 3", 4", 5°, 6" substitutions (52) and (53) 
into (51) obtained is system of equations for pressure: 

n+l n+l 



t F „ _ Pi+ij-Pu- l-sign(i-N x +\.5) 
n .M 2 

Xl+l 

1+ sign(i-N x +1.5) ! p ff -jw V 

2 ^-^-V,— T— ) X" + 



1 + ^(7-^+1.5) n+1 

+ H ^— < + 

~ t ~{ r vij+l T n+l 7 > ~ 






i =!,...,#, -1,7 =!,...,#, -1; 



.n+l „«+l 



— i(*uij- T n+i ; ) -<PuOjJ + (^7~)o / = 0, 

^i &ci ay 

1=0, y=l,...,iV y -l; 

ps i n+l n+l 

(° ( Pu\n+\ . \ UI7 n -T- Pn_~PiO\ ~«+i\_n 

C^ - )/o + T~ ' ^v/i _ r n+i : ) - (PviO i = °. 

y=0,i=WV x -l 
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3 rd scheme with central difference derivatives in continuity 
equation in nodes Q n 

n+1 n+1 n+1 n+1 

U i+lj ~ U i-lj V ij+l ~ V ij-l _ n ,, v ,, n/vw , 

The following schemes are used in boundary nodes 

t/" +1 t/" +1 n+1 n+1 , «+l , "+1 ,.«+! ,,"+1 

U Nj-^ Nx . lj V Nj -V N ._j «!■ -M . V . +1 -V . 

^ + ^1^ — ^- = 0, -^- ^ + ^- *- = 0, (55) 

» M+1 » n+1 v n+1 v" +1 ,," +1 _,," +1 v n+l -v n+1 

u iN y ~ u i-\N y v iN y - v iN y -1 "i+1Q M ;Q v /l v /0 _ n 

1 — 0, , z, 

by principle of interconsi stent approximation III central differences 
are used for pressure 

n+1 „n+l 



ij uij n+1 , , 

k ; ,,+h 



uij *„+l ^ , », (56) 

xj+I "•" r f*i 



n+1 n+1 

-,,"+! _ I?" ^- ^y +1 P'j~ l 

V ij -t'vij-^n^— —7— (57) 

Equations for p^. +l are derived by similar to [1] method: 

UF n _ T Pi + 2j-P« I- sign(i-N x +1.5) 

U uMj n+l h., 2 +hJ 2 

xz+2 xi+1 

l + «gw(i-iV r +1.5) „ +1 l + s/g»(/-1.5) _ P^-P"-L 

2 N ' j 2 [ui - l} " +1 /^ M j " 

1-^-1.5) B+ll 1 f l + «gnQ--iV y +1.5) n+ 

— ^n; i "t" I V w + 

2 ' h+h. + , 2 v 

+ ^+1 ^n+1 J 

n yj +2^ n yj+l ^ 



188 



n+1 n+1 ^ . / - i r-\ 

_ (J7 n _ Pij ~ Ptj-2 . 1 + sign Q - 1 . 5) 

K^ij-! T n+l ) 

n yj ^ n yj-i z ' 

l-sign(J-L5) +1 1 

1 v ,o h , — = ' i =i-,N x -l,j=l,...,N -1; 

boundary conditions for pressure are derived from (55): 

H xN x n xN x ^ H xN x -\ U J 

i = N x ,j=l,...,N y -l; 

hJ iF «> r " +1 h x2+ h xl } **' } + V°' °' 
i=0, j=l A^,-l; 

(-T-) «, + j— € Ny -( p viNy -i - ^ ~ h -^ h ) } - o. 

VX n yN y yN y T 'SW V -1 



J 



= N v ,i = l,...,N x -l, 



7=0, *=1,...,A^-1 



§5. Iteration algorithm. Global iterations method 

Solution of system of linear algebraic equations (48),(49),(50) in 

n+1 

relation to /?■ may be obtained by modified Jacobi method (simple 
iteration). For simplification of representing the algorithm upper 
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IS 



index n + l in p" + is omitted: p i . = p'? + . Number of ^-iteration i 
denoted as p f . For convergent iterative process 

lim^ =p?\v<uj i (58) 



Zero iteration p.. when k = on each time layer t n+1 is equal to 
value of pressure on preceding time layer t n : 

P°=p;, VtiJ I (59) 

k+l k ~ k ~ k+l 



pr- pj , i j F n - Pi*Th 

e k **" n+l h 



JI+1 
k+l _* 



1+ signji -1.5) _ P i} -Pi-ij . 

n, v r u i-lj T n+l , / 

2 »»■ 

l-5ign0"-1.5) „ +1 ■ i p*,-p* +1 ^ 

2 ' h iJ " +1 A ^ 

1 + sign(j -1.5) _ Py -Pij-i 

V vjj'-l 7 n+l , ' 



2* y*/— 1 rc-t-1 7 



l-jijB(;'-1.5) -iI+1 



v,; +1 =0, l</<^-l,l<;<^,-l (60) 



Iteration algorithm for boundary condition (49): 

k+\ k n+l n+l -. 

Pttf y ~PiN y U iN y ~ U i-lN y , J r n+l /7-n 

-^ + ^— + V '"' " ( '""-" 

fe+1 _ fe 
PiN y Pof-1 
~ T n+l —, )) = ( J, \<i<N t -l , (61) 

n yN y 
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n k+l A" — ^ n+l 

where coefficient with FiN . is equal to /i m ' 



KnM 



. Iteration 



algorithm for boundary condition (48) 

R+l 



k+l _ k k+l _ k 

PnJ PnJ 1 , „+l , rn PNJ Pn-1}., 



n+1 n+1 

+^ — ~ =°> l<J<iV -1, 

K 



(62) 



fe+i 4" _ »+i 

where coefficient with /? w ■ is equal to wj , , 

Here, < 6 < 1 - iteration parameter. In (60) collected are 
coefficients with pf in the form 



a;=t„j 



i 



1 I + sign ^-1.5 



-V 



xi V xi+l 



2h 



+ 



+ 



h. 



for explicit definition 
k+l 1 



1 l_+«gn4f-1.5 



s 



vVi 



2/r. 



6> /* . 



PLj 






(63) 



\ + sign(i -1.5) Plu \ -sign(i -1.5) iK+l -, 



1 V 

, IV- 1 vij 'n+1 , 



K 2 

l + «gnf-1.5 



M 0j J 

(64) 
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ij-1- 



(F^ l+ r n+l ^) 



1- sign(j -1.5) „ 



v I 



Similarly, from (61) and (62) computed are 



P iN y 






PiN u iN -tf-lNy —l v n+l _iF n 

- h 1 (TV,, xr viN v -l 



6 



+ T„ 



/?,. 



k 
Pff.-i 



l yN 



•y 



—) 



y 



p<i<N x -\ 



l yN y 



p k+l = 1 



e +*k, 



e h * NJ 



k «+l n+\ 

-(F n -1-7- P N *~ l J \ \_ Vn *J ~ Vn < - 1 i 



"}' 



(65) 



(66) 



h 



■xN r 



K 



>, \<j<N-\ 



Iterations (64),(65),(66) are continued until realization of continuity 
equation (24) with accuracy S « . Obviously, equations (24) and 
(50) are equivalent to each other because of representations (39) and 
(40), i.e. equation (50) is also a continuity equation, but represented 
in other manner. This circumstance allows to considerably simplify 

procedure of iteration algorithm for computation of /?■• . With this 

purpose (60) with the use of coefficient A-. is re-represented in the 
form 



k+\ k I k k 

(l + 0A") Pii ~ Pij + ( —4f"- t E^iZh 



M+l 
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I + sign {,-1.5^ 

™ ui-lj 



.ft „* 



Pa ~ Pi-u -> 1- signji -1.5) „ +1 

-*•-—£— - — 2 — "°' )+ 



T , X- 1 vij Si+1 , / 



"^ ' '"' / - 1 " +l h ( 67 ) 



v; 



l-«gnQ'-1.5) vB+1 \ =0 



Global iterations method for (24),(39),(40) in the first scheme 



Obviously expression in the left part (67), in brackets 

matches with (50), and (50) matches with continuity equation (24). 
Basing on this, it is convenient to introduce grid functions 



(68) 
(69) 



<■ 


, V tj by analogy with (39),(40): 




n+l* 

»,j = 


F n 

uij 


k k 
, Pi + U-Pv 

Km 




v," 1 ' = 


F". 


, pU~p\ 

Km 








at this 


ufyf 


are 


/^-iterations, 



accordingly, uf; -> v tj , that is, in the limit 
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lim «J +1 * = u£ +1 , lim vf k = v£ +1 , 
k — >oo /c — >oo 

since in iterative convergence (64), (65), (66) 

-, . k n+l \_i J • • ~~ 

Due to above-indicated match of expression in the left part (67), 
standing in brackets \ • • • / with (50), and (50) with continuity 



equation (24), contained \ ... /in (67) is replaced for iterative 
equation 

p™_p* u n+lk -u n+ } k 

' e k 



+ ^ V ^ ± - = 0,1 <i<N x ,l<j< N v 



(70) 



System (70) is equivalent to (67) and boundary conditions (65) 

and (66), is solved in relation to /?,^ + l , where by boundary 

conditions (26),(27),(28),(29): 
b" 



_1 nj / T «+l* _0 n+l* «+l* 

X — 1 + LI I Li, U n ■ — Z,U N _j . — U N _ 2 ; » 



n+l ^ n+l" n+l 

V 



^ = 2v^ _ ly . - v n Nx _ 2j ,0<j<N y ; 



on plate y = 0, w " +1 = 0, v," +1 =0,mk<i< mkk, (7 1 ) 

before the beginning and after the end of plate 

0n+\ k /a n+l* n+l* wo n+l* r\ 

."iO =( 4 ««1 -«i2 ) /3 > V «0 =0, 

1 < i < mk-l,mkk + l< i < N r , 
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y = H/L: w " +1 * = (4k "^ - m! +1 * 2 ) / 3 



v" +1 = (4v -v )/3i = 0N 

> > y 

Iteration process (68) - (71-73) is absolutely identical (64), 

(65), (66) and is realized all together until such value of k = k , with 
which continuity equation (24) in grid nodes is performed with the 

prescribed accuracy £ « : 

n+i k n+l k +l k +l k 

I u ij u i-lj v ij v ij-\ I „ 

| 7 + 7 \<E,V(i,j) (72) 

Iteration process is realized until criterion (72) is fulfilled. In all 
cases in satisfying inequalities (72) for solution on time layer t n+ j 
taken is 

k+\ _ n+l n+\ k n+ \ n+ \ k n+ l g .-> 

Pij -Pij ' U ij = Wy ,Vy =Vy ,Vij^ 

Global iterations method for (54)-(55)-(56)-(57) 

n+l 
when Py , coefficient in the scheme with central differences: 

A - T f 1 r 1 l-sign(i-N x +l.5) 

^, + i+^ hi + 2+h M 2 

1 1 + sign(i — 1.5) , 
+ ^^ -] + 

1 1 l-sign(j-N +1.5) 

+ [ - + 

Km+Kj h yJ + 2+h yj+1 2 

1 l + ^(7-1.5) ]} . = 1 ;! . = l ; 
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A -t 1 

KnSKn x +Kn x -i) 

A -r 1 

h yN r ( h yN y + h yN y -l) 

Ao =r »+i ,z' = l,...,iV -1 

+ h Ah 7 +h ,) 

Global iteration process is built as follows: 

ft * 

n+l* 77« „_ Pi+lj—Pi-lj 



ij uij n+l , , ' 

n* -n* 

n+l* = pn _ Pij+l Pij-1 

tJ ^ ^ h +h ' 

n yj+\ ^ H yj 

i = l,...,N x -l,j=l,...,N y -l; 

(1 + A.. 9) — + — + — — = 0, 

Ki+l+Ki h yj + l+ h yj 

l<i<N x -l,l<j<N y -l 

k+l k "+1 "+1* 

(l + A,/)^^ + ^-^ + (^)-=0, 

e h xNx dy 

k+l _ k n+l k _ n+l - 

(UAgjd) E^L + U JL^L + M K =0 , 



l<j<N y -1, 
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(i + A^)^-^ +( ^; + ^i=o, 

e dx h yNy 

k+l _ k p, n+l,k _ «+l 

0+iW) £fi_z£» +(%»' + "" "*» =o, 

# OK «,, 

Iteration process is topped when condition (72) is fulfilled 

u n+X u n+l v n+l v n+l 
I u ij ~ U i-lj V ij ~ V ij~l I „ 

1 h ■ h ■ 

n xi n yj 

n+\,k n+\,k n+\,k 
After that, last approximations V\] » "•« » "» are taken as a 

A* n+1 n+1 n+1 n+1 n+1 

solution P 9 *Pj- > M // * M i/ > V // *^ . 

Iteration parameter > may take any value from interval 
O<0<l/2. 

Global iterations method for (51),(52),(53) in the 2 nd scheme 

k _ k 



k _ k 
n+\ k _ r« _ fjjf ^J/'-l 

y vi)' n+1 



Jfc+1 _ fc n+1* _ n+1* n+1* _ n+\ k 
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0<i<N x -l,0<j<N y -l 

Global iterations for (54)-(55)-(56)-(57) by Krasnoselsky - Krein 
method of minimal residuals (m.m.r.) 

This is a universal method. Solution of a system of linear algebraic 
equations 

JV 

Y j a ij x j =b i ,i = l,...,N, 

7=1 

is found by iteration algorithm 

xf +1 =xf +O k R.,i = \,...,N , k = 0,l,...,k\ 
where R { - residual vector of /^-iteration 

N 

R i=H a , x j -b t ,i = l,...,N, 

7=1 

^=const - iteration parameter, is selected from condition of 

n 

minimization of residual k+1 — iteration V7i?* +1 ) 2 and is equal 

N N IN N 

k =-E£>,iW (La 9 R))*(^a,R)) 



'=1 7=1 / 7=1 7=1 



R K 



<S is fulfilled, 



Iterations are stopped when criterion max. 

l<i<N 

where < s - sufficiently small number. 

According to (m.m.r.) by (54)-(55) calculated are residuals of 
/^-iterations 

k k 

j.n+l" _ F n _ Pi+lj Pi-lj 

U ij r uij T n+1 j j ' 

n k - n k 

n+\ k =p n _ Pij+l Pij-l 

lJ Vtj ^ h +h ' 

i = l,...,N x -l,j=l,...,N y -l; 
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Rk _ Uj+ij Uj-ij v.. +1 v v _ t 



K + i+Kt h yj + i +h yj 

\<i<N x -U<j<N y -\, 

h xNx dy J 

K, = UjI ^H^<J<N y -i, 

n+l _.n+Uk 
nk _ ,°<P U \n+l . ™j V' k ,VH 1 

ox /z } ,j 

Thanks to specific character of the problem computed are 
accordingly 

J? k _ J? k 

nk . .__ '+1J '-!-/ 

nk __ n* 

i = l,...,N x -l,j=l,...,N y -l; 

after which defined are 

nk __ t-,£ t-,£ _ E>£ 



^+1+^ hyj + l+hyj 

l<i<N x -l,l<j<N-l, 
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AR k NJ =^f^,l<j<N y -l,A^=^,l<j<N y -l, 

-R k k 

AR^ y =—^,l<i<N x -l, AR* =^,l<i<N x -l 
V, h yl 

Iteration parameter is calculated 

N N y In N y 

(=0 .7=0 / (=0 .7=0 

and further pf = p\ +0 k tf,i = 0,...,N X , j =0,...,N y 
Iterations are stopped when inequalities are fulfilled 



l&*l= 



Kit -<i)/(h xi+l +h xi ) + (v;t -v£ )/(h yj+l+ h yj ) 
l<i<N x -l,l<j<N y -l 



<s, 



§6. Program on TURBO-PASCAL 

Major stages of programming the scheme of paragraph [1]. 
1 . Array declaration: 

i;^unIj^; +1 ^uniIjI 

v;=vnIj2>; +1 =vniIj: _ 

2*. Setting grid parameters N„ N y , h xi , i=l,..., N x ,h yj ,j=l,...,N y , T n+U 
included into equations constants Re,k, b° ,L,H. For Re»l it is 



U 

p 

p, 
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necessary to apply a nonuniform grid: h xl =h x2 =l/Re + , h x3 -2/Re + , 
h xNx =h xNx .!=l/Re + , h xNx . 2 =2/Re + ,h xi =(l-( h xl +h x2 +h x3 + h xNx + 
+h xNx .i+h xNx . 2 ))/(N x -6), 4<i<N x -3; h yl =h y2 =l/Re + , h y3 =2/ Re + , 
h yNt =h yN y. 1 =l/Re + , h yNy _ 2 = 2/ Re + ,h yj =( H/L-( h y ,+ h y2 +h y3 + 
+h yNy +h yNy . 1 +h yNy . 2 ))/(Ny-6),4 < j <N y -3, Re + = /VRe ,0 < y < 1 ; 
when Re<=500 analytical grid h xi =l/N„ i=l,...,N x , h yj =(H/L)/N y . 

3 . Setting of initial conditions (25) in cycle by i,j. 

4 . Setting in external (global) time cycle 



t n+1 = / t 1 ' m+ i of boundary conditions (26),(27),(28),(29) anc [ 

m = 

computation in cycle by i,j Q uij , Q vij , F uij , F V: 



m = 

in 
-uy » -^i-vy » ~ uij » ~ vij ■ 

5 . Programming of iteration algorithm 
(71),(72),(73),(74),(75),(76),(77). 

6 . Sending of calculated fields 

k+l _ b+1 n+1* „+l R +l* «+lw/^. 
Pg -Pij ' u ij = u ij > v ij = v ij ,Vij^, into arrays 

for u;=uNgjZ 9 v;=vNljZpi=PNtj-, 

Global iteration program for calculation of longitudinal flow past a 
plate with asperity by vicious incompressible fluid when 

H/L=8/JRs 
PROGRAM plastina; 

USES crt; CONST NX=140; NY=40;nxl=40;nx2=50;nyl=10; 
k=0;RE=40000 ; {in simulating turbulent conditions k#0, selected 
by experimentation } 

LABEL met, metk,met2,met4,met6,mettt,METPSI,MMM,SSS,TTT ; 
VAR outl , out2, out3 : text; 

UN,UN1,UK,VN,VN1,VK,PK,PK1,PN,PN1,AN,FUN, 
FVN,PSI: array[0..NX,0..NY] of real; 
HX: array [0..NX] of real; HY: array[0..NY] of real; 
KJ, i, j, KK, n, MK: integer; NN: real; 
HK, AM,AX, EPSILON, T, QNUij, QNVij, TAY, S, TET, 
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BO, HL, AMAX, UMAX, VMAX, PMAX, AU, AV, konu, konv, 
unn, aunn, vnn, avnn: real; 
begin 

assign(outl,'ainurl.pas'); assign(out2,'ainur2.pas'); 
rewrite(outl); rewrite(out2); 

TET:=0.01; B0:=0.1; HL:=8/SQRT(Re); EPSILON:=0.05; 
NN:=1;MK:=5; 
Writeln('hx[i] '); 

for i:=lto NX do begin HX[i]:=l./(NX-2*MK); Writeln('i= ', i, 
HX[i]:9:4); 
end; 

Writeln('hy[j] '); 

for j:=l to NY do begin HY[j]:=HL/NY; Writeln('j= ', j, HY[j]:9:4); 
end; 

{2 comment: setting of initial conditions. In simulation of turbulent 
conditions k#0 random disturbances are included into initial 
conditions } 

for i:=0 to NX do for j:=0 to NY do begin 
UN[i,j]:=k*sin((i+j)*0.25)/10;VN[i,j]:=k*cos((i*j !|! 3.46)/1000; 
PN[i,j]:=0; 

UNl[i,j]:= UN[i,j]; VNl[i,j]:= VN[i,j]; PNl[i,j[:=0; UK[i,j]:=0; 
VK[i,j]:=0; PK[i,j[:=0; PKl[i,j]:=0; 
end; TAY:= 0.001; 

n:=0; t:=0; 

met: n:=n+l; t:=(n+l)*TAY; 

for j:=0 to NY do begin UNl[0,j]:=exp(-B0/t)+k* sin(j*0.25)/40; 
VNl[0,j]:=0; UK[0,j]:= UNl[0,j]; VK[0,j]:=0; VN[0,j]:=0end; 
for i:=0 to NX do begin VNl[i,0]:=0; VK[i,0]:=0;end; 
for i:=l to NX-1 do for j:=l to NY -ldo begin 
if(i>=nxl.and.i=<nx2.and.j=<nyl) goto MMM;unn:= UN[i,j]; 
vnn:=VN[i,j];aunn:=ABS(unn); avnn:=ABS(vnn); 
konu:= ( aunn+unn)/2*( unn- UN[i-l,j])/ HX[i]+(unn- 
aunn)/2*(UN[i+l,j]-unn)/HX[i+l] 

+(avnn+vnn)/2*(unn- UN[i,j-l])/ Hy[j]+(vnn-avnn)/2*( UN[i,j+l]- 
unn)/Hy[j+l]; 

QNUij:=l/RE*(2/(HX[i+l]+HX[i])*((UN[i+l,j]-unn)/HX[i+l]-( 
unn-UN[i-l,j])/HX[i])+ 
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2/(Hy[j+l]+Hy[j])*(( UN[i,j+l]-unn)/ Hy[j+l]-(unn- UN[i j-1])/ 

Hy[j])-k*(( VN[i,j+l]-vnn)/ Hy[j+l]-(vnn- VN[i,j-l])/ Hy[j])/ 

(0.5*(Hy[j+l]+Hy[j]))) - konu; 

konv:= ( aunn+unn)/2*( vnn- VN[i-l,j])/ HX[i]+(unn- 

aunn)/2*(VN[i+l,j]- vnn)/HX[i+l] 

+(avnn+vnn)/2*(vnn- VN[i,j-l])/ Hy|j]+(vnn-avnn)/2*( VN[i,j+l]- 

vnn)/Hy[j+l]; 

QNVij:=l/RE*(2/(HX[i+l]+HX[i])*((VN[i+l,j]-vnn)/HX[i+l]-( 

vnn-VN[i-l,j])/HX[i])+ 

2/(Hy[j+l]+Hy[j])*(( VN[i,j+l]-vnn)/Hy[j+l]-(vnn- VN[i,j-l])/ 

Hy[j])-k*((UN[i+l,j]- unn)/HX[i+l]-( unn- UN[i-l,j])/HX[i])/ 

(0.5*(HX[i+l]+HX[i])))-konv; 

FUN[i,j]:=unn+TAY*QNUij; 

FVN[i,j]:=vnn+TAY*QNVij;MMM:end; 

if n>l then goto mettt; 

for i:=l to NX-1 do for j:=l to NY -ldo begin 

AN[i,j]:=TAY/HX[i+l]*(l/HX[i+l]+(l-(i-NX+1.5)/ABS(i- 

NX+1.5))/(2*HX(i+l))+ 

TAY/HYQ+l]*(l/HY|j+l]+(l-(j-NY+1.5)/ABS(j- 

NY+1.5))/(2*HY(j+l));end; 

mettt: 

for i:=0 to NX do for j:=0 to NY do 

begin PK[i,j[:=PN[i,j[; end;KK:=0;metk: KK:=KK+1 

for j:=0 to NY do begin UK[0,j]:= UNl[0,j]; VK[0,j]:=0; end; 

for i:=l to NX-1 do for j:=l to NY -ldo begin 

if(i>=nxl.and.i=<nx2.and.j=<nyl) goto SSS;UK[i,j]:= FUN[i,j]- 

TAY*((PK[i,j]-PK[i-l,j]) 

/HX[i]-k*(PK[i,j]-PK[i,j-l])/HY[j]); VK[i,j]:= FVN[i,j]-TAY*(- 

k*(PK[i,j]-PK[i-l,j])/HX[i]+(PK[i,j]-PK[i,j-l])/HY[j]);SSS:end; 

AMAX:=0; 

for i:=0 to NX-1 do for j:=0 to NY -ldo begin 

if(i>nxl.and.i=<nx2.and.j=<nyl) goto TTT;S:=(UK[i+l,j]- 

UK[i,j])/HX[i+l]+ (VK[i,j+l]- VK[i,j])/HY[j+l]; PKl[i,j]:= PK[i,j]- 

TET*S/(l+TET*AN[i,j]); 

if ABS(S)>AMAX then 

AMAX:=ABS(S);TTT:end; 
for j:=l to NY do begin 
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UK[NX,j]:=2* UK[NX-l,j]- UK[NX-2,j]; 
VK[NX,j]:=2* VK[NX-l,j]- VK[NX-2,j];end; 
for j :=0 to NY do begin 
PKl[NX,j]:= 2*PKl[NX-l,j]- PKl[NX-2,j];end; 

for i:=l to NX do begin 
if i<MK thenUK[I,0]:=(4* UK[I,1]- UK[I,2])/3) ; 
if i>(NX-MK) then UK[I,0]:=(4* UK[I,1]- UK[I,2])/3) ; 
UK[I,NY]:=(4* UK[I,NY-1]- UK[I,NY-2])/3; 
VK[I,NY]:=(4* VK[I,NY-1]- VK[I,NY-2])/3;End; 
For i:=0 to NX do begin 
PK1[I,NY]:= PKl[I,NY-l];end; 
for i:=0 to NX do for j:=0 to NY do begin 
PK[i,j]:= PKlfij]; UNl[i,j]:= UK[i,j]; 
VNl[i,j]:=VK[i,j];end; 
{ 1 6 comment: CHECK of fulfillment of continuity equation 

AMAX=max I «£ -uf)lh xM +(v^ -vf )/h yj+l \< e } 

If AMAX>EPSILON then goto metk; 

UMAX:=0; PMAX:=0; 

for i:=0 to NX do for j:=0 to NY do begin AU:=ABS(UNl[i,j]- 

UN[i,j]);AV:=ABS(VNl[i,j]-VN[i,j]); 

UMAX:= UMAX+AU+AV; 

PN[i,j]:= PKl[i,j]; if ABS(PKl[i,j])>PMAX then PMAX := 

ABS(PKl[i,j]); UN[i,j]:= UNl[i,j];VN[i,j]:= VNl[i,j]; end; 

UMAX:= UMAX/(NX*NY); writeln (out2, 'number of time layer 

N=',N:9:0, ' AMAX=' , AMAX, 'KK=' , KK); 

if N<500 then goto met; if N=1500 then goto met6; 

met4: NN:=NN+1; {17 comment: Ptintout each 1000 time layers} 

if NNolOOO then goto met2; met6: writeln (outl, 'number of time 

layer N=',N:9:0, ' AMAX=' , 
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AMAX, 'KK=' , KK );writeln (out2, 'number of time layer 

N=',N:9:0, ' AMAX=' , AMAX, 'KK=' , KK); 

writeln (' N= ', N:9:0); 

writeln ( 'Residual=' , UMAX, ' AMAX=' , AMAX, 'KK=' , KK ); 

writeln (outl, 'UN='); writeln (out2, 'PN='); 

for KJ:=0 to NY do begin j:=NY-KJ; write (outl, ']'=' , j); 

for i:=0 to NX do 

write (outl, ' ', UN[i,j]:9:4); writeln (outl); end; 
writeln (outl, '**********************'); writeln (outl, 'VN='); 

for KJ:=0 to NY do begin j:=NY-KJ; write (outl, ']'=' , j); 
for i:=0 to NX do 

write (outl, ' ', VN[i,j]:9:4);writeln (outl); end; 

writeln (out2, 'PN='); 

for KJ:=0 to NY do begin j:=NY-KJ; write (out2, 'j=' , j); 

for i:=0 to NX do begin PKl[i,j]:= PKl[i,j]/PMAX; 

write (out2, ' ', PK1 [i,j]:9:4);writeln (out2); 

end;end; met2: if NN=1000 then NN:=0; 

if UMAX>0. 01 then goto met; close (outl); close (out2); end. 
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Re=40 000, 140x40, t=tau*1 11 70 
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Field of averaged velocity vector V in longitudinal flow past plate. 
Re=40 000, 140x40, t=tau*1 0070 




Fields of pulsating velocity vector V ' in longitudinal flow past plate at 
various moments of time. 
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§ 7. Approximation of convection members without "schematic 

viscosity" 

Universal is multipoint approximation of a convection member, not 
containing "schematic viscosity", (ref. Jakupov K.B. 131), idea of 
which is as follows. 

From expansion into Taylor series 

3® h xt d 2 ® 

results 

3® ®" -®", h S 2 ® 
(—)» = lJ '~ lj + 2*L (^-) n +0(h 2 ) 

Standing here 2 nd derivative is approximated as follows: in near- 
boundary node 1=1 by formula 

tf<b ®" -®", ®", -®" ? 

cbc 2 ' ^ ft ft /* +/* +(AW * i+2+ ' W ' 

^ "xi+2 "h+1 '^i+2 + "ri+1 

In the rest of nodes i = 2,. . .,N x — 1 by another formula 
r) 2 ® ®" -®", ®",.-®" 2 

These approximations are applied in case of nonnegative coefficient 

u ;>o. 

From similar expansion into Taylor series 
derived is 

( & )j= ^ r ( &^« +0(A * |) 

Standing here 2nd derivative is approximated as follows: in near- 
boundary nodes z = 1, . . . , N x — 2 by formula 
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d 2 <5> 0" -O" O", -O" 2 

> u 2_ \« _ f '+2,/ ^i+lj _ ^<+lj 3j; x f_ 

fk 2 h h h +h 



( Vr) ; = (-^ — - - -^ — L ) . , + 0(^, +2 + ^ ) , 

l xi+2 ^xi+l ™xi+2 + ^xi'+l 

in near-boundary node i = N x — 1 by another formula 

These approximations are applied in case of nonpositive coefficient 
u" < 0. Discrete function sz'gn: 

for A<0, signA = -1, for A = 0, s/gnA = 0, /or ^ > 0, s/gnA = +1 

allows writing approximations not containing "schematic viscosity" 
in the form of a uniform formula 

ao \u"\+ u ; o"-o" 

,^ ^-^-iy ^-^ 2j ) 2 1 + ^/1(1-1.5) | 

2 ^ Vi Ki+Ki-i 2 

+ ( -) ]} + 

h h h +h 2 

n xi+2 n xi+\ n xi+2 ^ n xi+\ z ' 

n i n i -**« >f\n 



^ [( gW - Ql«j Q' MJ -_V[ 2 1 ~ signji - AT, + 1.5) | 

2 h h h +h 2 

^ n xi+2 " «+l "xi+2 T "xi+1 ^ 

+( ^ /, } ^ 2 H+00U- 

"« "ri-1 n xi rn xi-\ Z ' 

I « I , « Mini 

I W, ; I +«,, M.,-1 W„ I 



2 -<*£+ - 2 - q>^+Q(»;) (73) 

On analytical grid h xM = h xi = h x approximation is simplified: 
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„ = 3Qg-44>^ +^-2^1 + ^/1(1-1^) 
2h x 2 

3<t>; -20>Uj -W M j +®l2j A- sign(i -1,5) 



+ 



2h 



+ 



} 



o „ = ^j-Wt-^j^-signd-N, +U) | 
2/z x 2 

20; +1 , +2d>Uj -30; -0;_ 2J ^ l + signji-N, +1,5) 

2ft, 2 

Similar approximation of 2 nd order of accuracy is applied also for 

d® 

other convection members of v"( )" type. 

u dy 

Conditions of convergence and stability of explicit scheme impose 
constraints on time steps as follows 



d-Vi( 



3i«:i 



3 1 v" I 



Re/* 



■ + - 



2ft. 



■ + - 



■ + - 



Reft, 2 2h 



-)]>o, 



i = l,...,iV x -l,j=l,...,iV, -1, n = 0,1,2,. ..,JV r -l, 



/z x = min h . , ft = min ft , . 

i ■ J 

In algorithms of realizations of the above-given semi-implicit 
schemes in case of approximations of convection members without 

«schematic viscosity» changed will be only content Q\. , Q" r , 

included into F H ". , F" r . They will be represented with account to 
short representation in the form: 



Q" ■■=— I 



'^-K 



Re K M +h x 



Uy -U^ 



ft 



xi+l 



h. 



+ - 



K^K 
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( n n n n 



V »'+i 



3 



w y 



l«"l+«" <-l<l lv"l+v" v"-lv"l 

-{— ? L «"„ +^ ^u\ +^ y -u\ +^ ^u\], 

2 2 2 - v 2 



1 



fi" ■■=— \ 



Re ftrfrt+ft* 



^ v »._ v : v "-v" ,. N 



+ - 



***+** 



/ » n n n \ 



J 



\U- \+U- u - \U~ \V~ +v« v.- v» 

-{- ~v\ +~ -v\ + ^ *-v\ + ^ ^v"„} 

1 x i x i y i y 




Established field of velocity vector in cavity with the top moving cover. Grid 
140x80, Re=500. 
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§ 8. Semi-implicit schemes for solving of three-dimensional 
equations of vicious incompressible fluid 

For three-dimensional equations, suggested in 111 for simulating of 
both, laminar (k=0) as well as turbulent flows (k^O) of 
incompressible fluid 

du dp ^ d . du. d . du. d , du. 
P — + — = pF x + — (ju — ) + — (ju — ) + — (ju — )- 

dt dx dx dx dy dy dz dz 

d gl 2^ 8v 3^2 2^ 8w 

-k — [u% +w J(p-/Li—)]-k—[u\ +v J(p-ju—)], 

dy dy dz dz 



dv dp „ 8 . dv. 8 , 8v \ 8 . 8v^ 

> — + — = pF + — (ju — ) + — (ju — ) + — (ju — ) 

dt 8y 8x dx 8y dy dz dz 



; d r /2 , 2 >L du d r / 2 , 2>L dw 

-k — [v\ +w J(p-Li—)\-k — \y\ +v J(p-ju—)], 

dx ' dx dz ' dz 

dw dp „ d . dw. d , dw. d , dw. 
P — + — = pF z + — (ju — ) + — (ju — ) + — (ju — )- 

dt dz dx dx dy dy dz dz 

i d r y 2 2^*-, du ^ , d r g 2 2^, dv.^ 

-k — [wi 2 +w 2 J(p-jU—)]-k — [wi 2 +w 2 J(p-jU—)], 

dx " dx dy dy 

du dv dw 
— + — + — =0, 

dx dy dz 

d d d d d . , ,. . . 

— = \-u hv \-w — , with limiting conditions 

dt dt dx dy dz 

u \ t =o= d u (■*> y> z),v\ t =o = d v (*> y>z)> H t =o = <p w (*> y> z), 

u \cr=<Pu( x >y'Z' t),v\ cr=<P v (x, y, z, t), m\ a = <p w (x, y, z, t) 
Schemes and methods of global iterations given above are easily 
generalized. By analogy with directions x,y let's introduce nodes on z 



axis: 



z m ,m = OX..M z ,h m ,=z m -Z m _ l ,m = l,...,N z ,f(x l ,y j ,z m ,t n ) = f i ^, 
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Q h ={x { ,i =0,1,.. .,N x ,yj,j =0,1,.. .,N y ,z m ,m = 0,l,..N z }, 
Q h = {x i ,i = l,...,N x -l,y j ,j = l,...,N y -l,z m ,m = l,..N z -l} 

Let convection members in the general case be approximated by 

universal formula of (73) type, i.e. without "schematic viscosity" and 

let's take a two-dimensional scheme (57)-(60) as a sample. Let's 

n+\ 
introduce short notation of coefficient when Pij m in scheme with 

central differences for pressure gradient and continuity equation: 
A.. =r [ —*— [ — l -^n(i-N x +1.5) | 

m " +1 K M +Ki Ki+i+Ki+i 2 

1 1 + sign{i —1.5) 
h +h , 2 

xi XI— 1 

1 l l-sign(j-N +1.5) 

+ [ : + 

h yj+1 +h yj h yj+2 +h yj+1 2 

1 1 + sign(j -1.5) 1 

+ ^-^ -] + 

h +h . , 2 

1 r 1 l-sign(m-N +1.5) 

+ [ ^^ z - - + 

h m +h h x , + /* J _. 2 

zm+1 zm z/n+2 zm+1 

1 l + s7i?M(m-1.5) nl 
+ ]}, 

z=L...,A^-l,y = l,...,A^-Lm = L...,JV z -l; 

1 



NJm ^"n+1 



Kn(Kn+Kn-i) 
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A oj m =T n + i-rV 7 —J=l,...,N y -l,m = l,...,N z -l; 



A/V„m ^«+l 



1 



h yN r ( h yN y + h yN,-l) 



4-om =*-„ + i —: ; -— ,/ = l,...,^-l,m = l,...,iV z -l, 

h yl (h y2 +h yl ) 

A =r l 

iiN * " +1 h (h +h )' 

4,o=r„ +1 -— —,i = l,...,N x -l,j=l,...,N y -l, 

h zl (h z2 +h zl ) 

The global iteration process is built as follows: 

n+l s _r» _ t ri+ljm Pi-ljm K n / nl 

U ijm ~ * uijm ~ '? n+li ^ ^ ■" ~| ^ \- U ij+\m \ U ij+\m ~< 

h xM+Ki Vl+«,j 

. n2 \-0,5 s n / nl nl \-0,5 .v i . 

+ W jW Pij + \m- Uij-lmWij-lm +W iHm) Pij-\mi + 

_i_ / T w ( nl . nl \-0,5 5 

+ K l U ijm+l\ U ij m+ l + V y m +1 J Pj/m+l " 

n / nl nl \— 0,5 s ~l // / , / \1 

- M i/m-l ( M i/ m -l + V i/m-l ) P,jm-1 J '("zm+1 + \n ) J > (74) 



«+l S —J?" _ r f -PiZ+lm "i/-lm "- r n / nl 

ijm vijm n+ll » » » » L i+\jm \ i+ljm 

n2 \— 0,5 .v rc / nl nl \— 0,5 .v n 

+ W /+l> ) A+lym - V i-l./m t V /-Um + W i-1 jm ) A-i> J + 

i r n s nl nl \— 0,5 .« 

"*" ^ L V ym +1 V M i/m +1 """ *ym +1 / Pijm+\~ 

~ V ijm-1 (KjLl + V £-l )^' 5 Pjn-1 ] Kh m+1 +h 7m )}, (75) 



n+l s _ T-m r rijm+l rijm-l "■ r » A, "2 . 

"zm+1 + "zm ft xi+l + 'lei 
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+>C)^i>. - <^U + <^' 5 pI^ + 



, j r n / nl . nl \-U,3 .v 

- <_ lm {uf_ ln + wf_ lm )-°' 5 pl_ im ] /(/Z w . +1 + h yj ) } , (76) 

i = l,...,N x -l,j=l,...,N y -l,m = l,...,N z -l; 

p s+l -p s . u n f. -u n f. 
(] j. A R\ — — ' i — — m i 







h xM+ h » 



v.. , — V-'i W»", — W- . 

ii+lm ii-lm ipn+l iim-l ^ 

+ — — + — — = 0, 

i = \,...,N-\,j = \,...,N-\,m = \,...,N-\; 



S+l J 



n+1 n+1 



(1+4^) 



PNJm~PNJm ^V> "j^y, 9^ „ +1 %>. 







■ + 



n+1 



/? 






.vJV, 



dy 



d: 



'NJm 



s+l _ s ji n +' — rn n+ ^ a a 

/,,! ^ FOjm lojm M ljm V^Ojm .^Cff „ +1 C^> n+1 

U + A), ffl #) " + , + \—)opn + \^>Ojm = U > 

9 h rt dy dz 



"x\ 



\<i<N-\,\<m<N-\\ 



.s+l 



.n+1 n+1 



d+V.*) 



(i+Ao^) 



PiN y m PiN y m ^viN t m V i5V,-lm fity „ +1 (3<^ „ 



6> 



-+ 



A 



/-r^^+i + rr^\«+i _o 



.s+l 



ftp*-/ ? 
(9 



n+l s n+1 



ck 



& 



iOm , ilm 



■ + 



PviOm Jf., + 1 , A 



n+1 



A 



T V ~ /iOm T V ~ A'Om u ' 



vl 



& 



dz 



l<i<N-\,\<m<N-\; 
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S+l _ S /y)" +1 _ vu" 

(l + A iW _60 + + K-r-)m,_ +\.-r-)mi_ = u > 



6 h 7N dx l]Nz dy 



S+\ s , n+\ s rn «+l 



U + A/o 6 ') — " — + ; + ^^'o + ^~></o = U ' 

6 h zl dx dy 

l<i<N x -l,l<j<N y -l 



Global iterations of finding pressure by Krasnoselsky - Krein 
method of minimal residuals (m.m.r.) 

By equations (74) - (75) - (76) computed are three components 
u yl » v lpn > w yl > a f ter which, residuals of ^-iteration 

„+!« n +l S n+l s n+l s n+\ s n+l s 

jps _ i+ijm i-ljm ij+lm y-lm ijm+\ ijm-l 

' Jm K + l +Ki h yj + l + h yj K + l +h w 

i = l,...,N x -l,j=l,...,N-l,m = l,...,N-l; 



_ < PuNJm U Nx -ljm ,d<P v y+l /^ z \n+l 

nj _ M ljm ~ ( PuOjm , 0<P v . n+1 , 0<P z .„ +l 

K 0im ~ , + V _ /0/m" 1 "^ ~ /0/m' 

^,i dy & 

1 < / < N v - 1,1 < m < N - 1; 

n+\ _ n+\ S 
, _ ( PviN y m V iN y -lm 0$ B+1 ( Q(p w , n+ \ _ „ 
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1+1 o o 

ns _ V ilm ~ ( PviOm , /^«\«+l , / ^w \»+l _n 
A iOm - , """V - ZiOm ""V - /(Om _ U ' 

w yl ox <3z 

\<i<N -l,l<m<N -1; 

ns _ jW. "^V,-! /^nxn+1 /^yxn+1 

1 < i < N x -1,1 < 7 < iV y -1 

Thanks to specific nature of the problem computed are accordingly 



R s — R s 

p.' _ f i+ljm i-ljm , k 2 

A «ym --^+li— " — + —[U" j+lm (U" j+lm + 

KM+Ki h >jH+ h yj 

. , ,n2 \-0,5 D .s n / nl , ,.,«2 \-0,5 jV 1 , 

+ W i/ + lJ R ij + lm- U ij-lm\ U ij-lm + W ij-lm) K ij-lm\ + 

. a r n i nl . nl \-0,5 p.s 

"I 7 ^ W" + 1 ^ M ')' m + 1 IP" + 1 ' >J m + 1 ~ 

"an+1 """"an 



/" (w" 2 +v" 2 ")~°' 5 /? s 

i/w— 1 V ijm—1 ijm—\' ijm- 



p« __ r '>' +lm 'J- lm i r,," j-,," 2 , 

^vijm ~ T n+\\ , h h h *- V i+ljm\ V i+ljm ~*~ 

rc2 \— 0,5 T>s n / nl nl \— 0,5 7>v l 

+ W ;+lim ) ^+1 j» - V i-\jm V V i-\jm + ^-1 jm ) "i-ljm J + 

"7 7 L pn +1 V M ; )m +1 "•" V ijm +1 / ^j/m +1 ~ 
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v n . Au n2 ,+v n2 ,)^ 5 R S . .]}, 



r" 2 ,+v" 2 J" ' 5 /?.* 

D» __ r f R ijm + l- R ijm-l * F " / «2 

V" W, , +, , iyy Mjm \y Mjm ^ 

n zm+l +n zm n xi+\^ n xi 

. nl \ -0,5 r,.s ,.Ji A.B2 ,.,,12 \ -0,5 pi n . 

+ w 1+Um ) R Mjm - w iAjm (y iAjm + w iAjm ) K lAjm J + 

, "■ r n / nl nl \— 0,5 tjs 

+ ~ —^ W ij+lm\ U ij+lm +W ij+lm) K ij+lm- 

h yj+\ + "»■ 

n / nl . n2 \-0,5 r>,? n i 

- W ij-\m\ U ij-\m +W ij-lm) K ij-\m\i > 

i = 1,...,N X -1, j =l,...,N y -l,m = l,...,N z -1; 

after which defined are 



.v _ ui+ljm ui-ljm t vij+lm vij-lm t wijm+1 vijm-1 

1,1,1 &• , , + h : h ; , , + h ; k m: i+ h 



a j-j.v ui+lim ui-ljm vii+lm vii-lm 

AR i im = — : ; + — : ; + - 



xi+l xi yj+l yj zm+1 zm 

i = l,...,N x -l,j=l,...,N -l,m = l,...,N -1; 



>■ 



AK xJm =~- R f^ L ,ARl jm = ^,1 < j < N y -1,1 < m < N m -1; 



Kn x k 

h y\ V, 

A^ =^,A^=^^,1</<A^-I,l<j<^-1 



Iteration parameter is calculated 

(=0 ;=0 m=0 / i=0 ;=0 m=0 

andfurther^ 1 =^+^;,/ = 0,...,^,;=0,...,^^ = 0,...,^ 
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Iterations are stopped when inequalities are fulfilled 

\R; m \<s,i = 0,...,N x ,j = 0,...,N y ,m = 0,...,N z 

* * * 

s* n+l S n+\ S n+\ S 

Last approximations />- m , M-. ,V- ,W-- are taken as decisions 



«+l n+l 



n+1 



p.. » p.. ,«.. ««.. ,V- »V- ,W- »W- , 

r i;m r i;m ' " nm " am ' i/m wm ' i/m am ' 



i/m r ijm ' I/m 



<;m > ijm 



ijm ' i/m 



ijm 



difference approximations of expressions 

d , du w d , du w d , du w 
— (ju — ) + (jU — ) + — (ju — ), 

dx dx dy dy dz dz 



k — [u% 

dy 



+ W JjU — ], k — [w| 



2 2 ^l 8W 

. „ +v ^ 2 // — J, 

9y & dz 



d , dv. d , dv. d , dv. 
— (jU — ) + — (ju — ) + — (ju — ), 

dx dx dy dy dz dz 



ox ' ox dz ' oz 



d . dw. d , dw. d . dw. 
— (jU — ) + — (jU — ) + — (jU —), 

dx dx dy dy dz dz 



k — [w\ 

dx 



2^- du 1 , d r / 
+ w JjU — ],k — [w\ 



2 ~=4 ^V i 

dx dy dy 



are included accordingly into F^ jm ,F v " jm ,F" ijm , it is advisable to 
approximate them by samples (ref. 121): 



d_ 

dx 



A 



SO 

dx 



h +h 

_i ijm xi+l xi 



2 p y^i+ljm "■" ^ijm X^i+ljm ^ijm) 



2h 



'xi+l 



i v 1 -\- 2 n VcD" — db" 

^ ijm i—ljm ' ^ ijm i—ljm 

2/j. 



]+0 (h xM -h xi )+o(h 2 j, 
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d_ 
dy 



A 



dy 



u"i 



h yj + l +h yi 



w^+^mi 



O" ) 

ijm ' 



2h 



W+l 



a n _i_ ji n Vd)" — d)" "i 
ijm ij-lm '^ ijm ij-lm ' 



2h.. 



] + 0(h yj+l -h yj ) + 0(h 2 yj ), 



d_ 

dz 



A 



SO 

dz 



K + i +h 



-[■ 



(^"iim+i + A" m X^L+i - ^L ) 



(^■L +^L-i)( ( i ) L 



2/j 



zm+\ 



ijm — 1 



2/j 



l+W^-O+^O 



These approximations are suitable for discontinuous coefficients /5/. 
Afote. Generalization of stated iteration algorithms by three- 
dimensional problems in cylindrical, spherical and other systems of 
coordinates does not cause special difficulties and is performed 
similarly to the above-stated (ref. III). 
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Chapter 10. PARADOX OF COMPUTED "NEGATIVE" GAS 

DENSITY. 

EFFECTIVE METHODS OF SOLVING COMPRESSIBLE 

GAS EQUATIONS 

§1. Paradox of negative gas density is caused by incorrect 
realization of law of mass conservation 

In numerical solution of equations of compressible gas by explicit 
difference schemes of Brailovskaya, Davydov, by implicit schemes 
of III, 111, 141, 151, IY1I and others, obtained are negative values of 
gas density p , which conflicts physics of this value, measured in SI 
by relation of kilogram by cubic meter (mass of any matter is 
positive). Reason for this is incorrect realization of difference 
analogue in continuity equation 

— + divpv =0 (1) 

dt 

used in the above-mentioned works for direct calculation of density 
p by difference scheme of the following type 

rc+l n 

A^ ?- + (divp n v n ) h =0, (2) 

where A - linear invertible operator; (divp n v n ) h — difference 
approximation divpv . Indeed, from (2) calculated is density on the 
following time layer by inversion of operator A : 

P =P ~T n+l A (divpv ) h , (3) 

It is obvious, that in those points of flow, where A" (divp n v n ) h 

preserves in the counting process a positive sign, density p n+ with 
the growth n will be decreased and will become negative, which is 
absurd, since it is deprived of sense, or vice versa, increase up to 
unreal large values, if A" {divp"v") h >0. 

G.I.Marchuk, for eliminating appearing fluctuations increasing 
when t — > co by amplitude, even in observance of conditions of 
stability, proposed density p profile smoothing on each time layer 
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(first towards x, then by 0) by formula p = p m +(p m _ l -2p m +p m+l )S, 
and for stability of counting it is necessary that < 8 < 0,25. 

These facts and a number of other facts brought up question of 
fallacy of density calculation by algorithms of type (3) and search of 
in principle other algorithms of solving equations of viscous 
compressible gas. 

Such algorithm was grounded in /6/,/14/ and this book, where is 
was concluded that realization of law of mass conservation (1) is 
performed via pressure function 

8 2 . p , ^,8 2 p 8 2 1 

3t 2 RT tT dxf dxf 3^ ^ 

8 
+ — [div{pv i v)-div{jugradv i )-pF i ]}, 

ox i 

i.e. closed system of difference equations for calculation of pressure 
is constructed from difference analogues of continuity equation. 
Having in mind that pressure gradient is included into viscous fluid 
equation, both compressible as well as incompressible equally, only 
for compressible gas density is defined by pressure and gas 

P 

temperature p = , and for incompressible media p = const , 

technology of constructing these finite -difference equations is 
absolutely similar with the technique stated in this book for 
incompressible fluid equations. 

In monograph 161, with the use of method energy inequalities. For 
the first time studied were issues of stability and one-valued 
solvability of family of semi-implicit and semi-explicit schemes for 
equations of Navier-Stokes, thermal conduction and diffusions with 
account to chemical reaction of burning in cylindrical coordinates 
system. Closed system of difference equations was obtained for 
pressure. From foreign works on numerical solution of compressible 
gas equations, it is necessary to notice article of D.G.Lilley "Simple 
method of calculation velocities and pressure in strongly swirling 
flows" // Rocket engineering and astronautics. 1976, v. 14, June, p.57- 
67. This work is close in main idea but in it, in three-dimensional 



221 



problems used are 4 half-step displaced in relation to each other 
spatial grids, in two-dimensional - one grid less, while in the 
methods of the given chapter one grid is used, which considerably 
simplifies theoretical analysis and programming. Half-step displaced 
grids contain values of unknown functions beyond physical domain 
of flow, and their definition introduces uncorrectable errors in results 
of calculations, which is discussed in details in Chapter 8. 

§ 2. Semi-implicit differential scheme of complete equations of 
compressible viscid gas 

Principles of building of semi-implicit, semi-nonimplicit, method 
of rhythmic steps and variable directions of schemes for equations 
Navier-Stokes of perfect gas were expounded in monograph 161. 

In this paragraph technology 161 is applied for new system 
equations, well-founded in Chapter 1: 

r dV e , V 1 dV e\ , d P Z7 ■ 

P { ^T + L V >n ^ )+ ^~ = P F e + 

ot m=i ax m dx e 

+ f j ^( J u^-)-^[(\ J u- J u')divJ], *= 1,2,3, 
m=i dx m dx m dx e 3 



,dT J^ dT 



P°v ("T" + 2 v - 7T - ) = div(AgradT) + 
ot m= i ox m 

e=l m=l C^m J 

7T + S^ = 0, P = /JyriT,// = JU(T)A = MT), 

at t^i «* e 

under entry conditions 

r L=o =d T >f\t=o = d P ' v e\ t =o =d e ,e = l,2,3, 

and boundary conditions on boundary of district of integration 

T \ S =<Pr> V e\s=<Pe> e=1 > 2 > 3 
Notice, that for pressure p entry and boundary conditions don't put, 
for density of gas p don 't put boundary conditions, there are not 
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necessary. First, offer semi-implicit differential schemes are 
organized thus, that values of pressure are calculated in appropriate 




Phc.1 

boundary node S h of net domain by iterative numerical algorithm. 

Therefore density of gas is calculated in full field, including 
boundary, from equation status p — p/(RT) . By way of example on 
drawing 1 , list problem statement about convection by heating of gas 
below T\ s = (p T , rest walls of field are heat-proof. Clearly, that in 

this task any assignment of boundary conditions for pressure will 
appreciate absolutely other task. 

Secondly, for numerical solution of these equations completely 
apply technology of building and realization of schemes, stated in 
previous chapter 8. In technology of chapter 8 for pressure proved 
difference boundary conditions from continuity equation. The same 
circumstance takes place and for equations of compressed gas. 

In tasks of flow of bodies unlimited stream possible decide 

complete equations in net domain < i < N l ,0 < j < N 2 , < k < N 3 
put, for example, in tasks of calculations of flow by direction x , by 
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way of one from boundary conditions by i=N] equality of zero of 2 
derivative from desired function 

d 2 f/8xf=0 (1) 

On nonuniform mesh (1) approximate by formula 

[(/£ -O/ V -ift-^-f^'K.^'K + V- 1 )=a (2) 

on analytical grid /z 1A , = /z 1A , _j = fy from (2) prove 

NJk - Z 'Jn 1 -Ijk ~ JN l -2jk ' (3) 

Approximate all convection terms, including 

(vf3v; / dx x \ Jk ,e = 1,2,3, (y»&T» 1 3x^ ijk 
possible make by formula without "diagrammatical tenacity" 
„ .dm. \v n J+v n m , v n m -\v n m \ „ 10 „ 

V mijk (^")* = 2 ^ 2 ^ i " VmCO ' X ™ ' m = ' ' ' 

or use principles from monotonous schemes of chapter 6 in some 
cases by more simple scheme 



da. 



v„ 



+ v v -■ 

but in any case with central difference 

„ da „ o)- +m 

v mijk (—)ijk = v m — "^ — "S m = i- 2 - 3 , 

m 

they let to incorrect results (to serrated decision). 

In semi-implicit scheme gradients of pressure and component of 

speed in continuity equation approximate with central differences 

and take on upper layer n+1 outward iteration, but all rest terms of 

equation on lower layer of time t n : 

3 
-«r/_.n+l ,.ii\ /_ , X ' n n n , _n+l 



P L(v e - v J / r„ +1 + 2^ v m v &m ] + ft = p F e + (4) 

m=l 

3 3 

+ Z ( '"" V ^ ) * m "E ( '" n/3 "'"" ,)V ^ ) 5:J' e = 1 ' 2 ' 3 ' 
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p\[(T n+1 -T n )/r j+1 + Yv n J.l] = ^(AXX + 

m=2 m=\ 

+M n j]i(v:o 2 -p"i^k-(\M-M'")(j]v:~ x f (5) 

Schemes (4), (5) noticed for indexes 1</<JV, -1, \<j<N 2 -1, 

l<fc<iV 3 -l. In particular, possible use monotonous schemes. In 

an effort to simplification of statement applied approximation 
"without diagrammatical tenacity". 

1°. Scheme of continuity equation 

In schemes (4) gradients of pressure are replaced by central 
difference derivations. By stated in 161 principle coordinated 
approximation, derivations in continuity equation also approximate 
by central difference derivations'. 

« -£)/*„. +(A n+ % HpX% 2 HpX% =o, (6) 

i = l,..N v j = l,...,N 2 -l,k=l,...,N 3 -l 
In boundary nodes assigned temperature and speed components 

= q> n e +l ,m = 1,2,3 



rpn+l 



n+\ n+l 



h 



In boundary nodes i = 0, V(y, k);j = 0, V(z, k); k = 0, V(/, j) 
continuity equation approximate by technology III difference 
forward 

tfHpXX +i^^ = 0,xzS eh ,e = 1,2,3, (7) 






<9x„ 



in boundary node i = N v \/(j,k);j = N 2 ,\/(i,k); k = N 3 ,\/(i,j) 
k = N 3 ,V(i, j) approximate differences rearward: 

p:Hp n v n :\ +Z^=a«v =u,3, (73) 

For example, in boundary nodes j = 0,\/(i,k) (7) transform to type 
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(Pm ~ Pm ) l Vi + (Plmkfmok ~ PlmVtm ) UK* + K ) + (/Mi " (8) 

-A %" t) //l 2i +(/4 + i?fi + i - Pm-Aa) '(K* +K) =( W =o,v(U) 

In boundary node A: = 0, V(z, _/') notice (7) take type 

(A"o' - A"o ) / T n + 1 + (PMJ0<PIMJ0 ~ Pi-ljo9u-lj0 ) %tl +K) + (A" + 10<!l0 ~ 

-p;^<p^o)i(h m +h v )H P ;^-p; ^/h n =o^=oMj) 

In boundary nodes i = N 1 ,\/(j,k) notice (7a) take concrete form 

/ —"+1 „" \ / _ i / —« ™"+l ,," ,.n+l \ / z„ i 

(P^/fc -PNJkl'tn+l+iPNJkPlNJk -P Nl -ljk V lN 1 -ljk) / 'hN 1 + 

+ (PNJ+UcPwlj+lk ~ PNJ-lkPlNJ-lk ) /(^y+1 + Kj ) + 

+(p;ui<* + i-/4w<*-i)/(^ + i+*3*)=o,j=av(a) 



2°. Realization of state equation p - pRT 

By this clear scheme (5) count field of temperature T n+l , 
therefore state equation possible present in type p" = p" /(RT n ) , in 
difference continuity equation (6), (7), (7a), (8), (9) and other put 
Pm 1 = Pijk 1 l(,RT£ l ), because p" well-known value of density 
with previous layout of time t n . It necessary take into account, that 
state equation p — RpT take place in all points of net domain 
< i < N v <j< N 2 ,0 <k<N 3 . 

3°. Difference equation for pressure 

Previously equation (4) increase on step by time T n+l , after that 
present in type 

n n+\ r «+l ^-\« -i 

Pijk V eijk ~l~Px e ~ { ^eijki T n+V e =1*2,3, 

3 3 

eijk=Pijk\-- T n + l V eijk +L V myk V ex m ~ F e ^Z/^-A + 
m=\ m=l 
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+X (» n /3-»' n )vlJ- Xe , (10) 

m=\ 

l<i<N, -1,1 <j<N 2 -1,1 < fc < N 3 -1 

With purpose of using of results previous chapter come back to 
initial indications 

v x =u,v 2 =v,v^=w,(p x =(p u , 

<p 2 =<p v ,<p 3 =<p w ,h 1 =h x ,h 2 =h y ,h 3 =h z ,N l =N x ,N 2 =N y ,N 3 =N z , 

F" =-t O n F n =-t O n F n =-r O" 

In this indications (10) take type 

n+1 _ „+i 

n n+1 j^n ri+\jk ri-\jk 

Pijk U ijk —f uj j k —T n+l 



Km 


+Kt 


? 


n+1 
Pij+lk 


n 
'Pij 


+1 
-lfe 



n n+1 _ 77" _ J- »+ik J- i] -IK 

Pijk V ijk ~~ r vijk T n+l , , ' (11) 

n yj+i + n yJ 

n+1 n+1 

n n+1 77 n "ijk+1 *ijk-\ 



n" w = F n —T 

rijk ijk wijk n+1 » » 



Z&+1 ' ""zS 
n+1 

Equation for pressure /?» prove by way of substitution indifference 
analogues of continuity equation (6), (7), (7a), (8), (9) and other: 

[p;: l /(R^ l )-^ k yr n+l + 

4f „ _ Pltijk - PT , 1 - signji -N x +1.5) 

' ~\ r ui+ljk ^n+l , , > ~ ' 

h xi+ 2+K M 2 

n-\-\ H+l 

1+ signji -N x +l.S) „ +1 1 + signji -1.5) _ P ijk -Pi-2 jk 

"■" r, 'PuNJk r, y^ui-ijk ^n+1 , , > 

l-«g»(/-1.5) , +1 -, 1 , }+signjj-N y +1.5) 

2 ^ JZ+~K i 2 ^ 
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^ (Fn _ t P^ k -p;;\lsignU-N y+ l.5) 

+^ vij+ lk T n+ 1 h h ) 2 

n yj+2 ^ n yj+\ A 

n+1 n+1 ^ , . * »» x 

_ (J7 n _ Pijk ~ Pij-ik , 1 + signjj - 1 .5) 
l/Vu T »+ l h +h ) 2 

1- sign(j -1.5) n+i ^ 1 

2 — % '« ^-+^: + <i2) 



w'+i w 



J+s/gn(£-Af z +1.5) ! 

+ "T " ; <P wij N z + 

(J7 n _ PT+2 ~ Pijk' v 1 ~ g^g Z jV g _+l -5) 

h zk+2 +h yj+l 2 

n+1 n+1 -t . /•, -i /-^ 

_ (F n _ Pijk -Pijk-2 , 1 + sign (k - 1 .5) 
1 -sign(k -1.5) B+ i -i 1 _ n 

2 \k + l+Kk 

i = l,...,N x -l,j = l,...,N y -l,j = l,...,N y -l, 

boundary conditions for pressure prove in type: 



[pZk'(RT^ k )-p n NJk VT n+1 + 

. n+1 n n+ l 

+ -7—WuNJk ~ ( F uN,-Xjk ~ ^+1 —T 7T ) J + 

"■xN x ' l xN x + ^xA^-l 

+ (^ft + (^)& = 0J = U,N y -lk=l,.,N z -1; 

[p"ot W? ) - a% V *»*+-£-{ (Kij k - 
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n+1 „n+l 



' 7 n+l 



- 7 « + i -; ; — ) - (Puajk ) + ( — - — hjk + ( — : — )oj» = °> 

./ = l,..JV y -U=l,..JV z -l; 

[ Aw) /(^^l ) - /&,* ] / 7 „+i + 7 — %?,* -(C y -u - 

i = l,...^-U=l,...^-l;, 

h yl h y2 +h yl 

p. n ..n+l Flri n /r> n+ ^ 






-^,o}+( a )+( a ) = 0,i=l...,N x -l,j = l,...,N -1; 

§3. Global iterative method of calculation of pressure 

In this method in iterative process is included equations (10) on 
every iterative layout by s: 
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n s,n+\ r i,n+l s^n -i 



s,n+\ 
eijk 



S ~ ( Teijk ' 



(13) 



e = 1,2,3,1 < i < iV, -1,1 < ;' < N 2 -1,1 < Jfc < N 3 -1, 



which in complex type assist in difference continuity equation: 

(i+^)(^" +1 -pr 1 )^+[pr 1 /(^ +1 )-4]/T; +1 + 






ijk ' *-*ijk ^ ijk ' 'ijk - 



where easily calculate from (12) 

1- sign(i-N x + 1.5) 



(14) 



A". = ■ 



1 



■+[■ 



1 

x/+2 jci+1 



+ 



+ 



' 77+1 ijk 

1 l + sign(i-l.5) r 



Ki+Ki-i 



-1 



Km+K 



■ + 



+[■ 

+[ 



h yJ+2+K +1 

1 



l-«gw(j-Af y +1.5) 1 l + «gn(y-1.5) 1 



A «+Vi 



Ki^+K 



\-sign(k-N, +1.5) 



1 l + «'gn(£-1.5) 



K + 2+K + i 



K + i +h zk 



2 **+Vi 2 

1<7 <iV, -1,1 < j" < iV 2 — 1,1 < it < iV 3 — 1 
Iterations in boundary nodes j = : 

+ #W Xfto* " Aot )/ + [Pm !( RT m ) " Ao* 1 7 T « + i + 



+ (A+io*Cio* - Pi-m^-iok ) ^i.+i + ^ ) + (A> 



n .v.n+l 

' 7 2,u - (15) 



n+K 



-Pm<P2m) lh 2\ +(/4+iCow -Pm-A-i) f (h ik+l +h 3k ) = 0, 



Aok - 



T n + l RT i"ok 



- + - 



^21 ^22 "■" 21 



,/ = l,...,A^ 1 -l,^ = l,...,^3-l 



Iterations in boundary nodes k = : 

a+ayoC" 1 -PiT)io+\^ l i{Rj;; l )-^ir n+l + 

y ij+WY2ij+W '(16) 



+ ( A" i /oftw 70 - A"i ,o0u-I ,o ) AA W + *u ) + (pLoPm 1 



MjOYK+ljO 

rt+1 



-lj0Y\i-\j0J 



pLoC-; )% + i +w+(A>sr 1 -&<0'K =°> 
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A" 1 I r "+i 1 i -1 A/ -1 r-1 N -1 

l n+l 1 ^ 1 ijQ "31 "32^ "31 

Iterations in boundary nodes i = : 

+ (A^„^ -pOjMoj^/Ki+iPoj^PlOj+ik- ( 17 ) 

~ Po H k<P2 HkV( h 2] + l +h 2j) + (PQjk + l<Pm + l ~ AVl^oi-lV^i+l + M = > 

4" = 1 +^_!_,y=l,... J tf 2 -U=U,tf3- 1 

Iterations in boundary nodes i = N l : 

(i + m;_ ., x^r 1 - *o i o + [<;; +l krt^ ) - ^ ] / Vl + 



2 ' 



+ yPNJkYlNJk ~ PN^ljriN^ljk > ' ^Ny + ^PNJ+imNJ+lk ~ V i 5 -' 
-plj-lkVlNJ-lk) 1 ^! + } h}) + (p n N l m$N l jM -pljk-lfPiNJk-d^M + M = > 

l 'n+l lyl N 1 jk "liV, "liV, ^''l/V.-l 

Iterations in boundary nodes j = N 

a+^)wr-^^+wr / (^)-^ / vi+ 

+(#V<U* - pliN 2 k<p"^N 2 k) /(K + i +K)+(p'k_k ( Pm 2 k - (19) 

-Pw 2 -u v 2^«)/^ 2 +(P^ 2k+i ^ lM -Pm 1 kA^ 1 k-i) / ( } hk + \ +K) = > 
A" = — r + ^i , l</c<AL-l,l<z'<Ar,-l 

tfjt ^n+1 , , , ' 3 1 

l, n+l lyl iN 2 k n 2N 2 n 2N 2 T ,l 2/V 2 -l 



Iterations in boundary nodes k = N 

-A"-l« 3 ^- 1 l W3 ) / ( /j 2; + l+ /j 2;) + (^3^3 ~ P^-l^-l) 7 V = > 
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Aw = IT + — ,\<i<N,-l,l<j<N 2 -\ 

' l3 x RT" + h h +h 

Parameter of iteration select from interval < < 0,5 . In all 
without execution applied above iteration processes, exactness 
calculation of field of pressure limited by exactness of execution of 
continuity equation. 

§4. Indexed family semi-implicit schemes 

Set above scheme is special case of family semi-implicit schemes, 
well-founded in 161. 

3 _ 

n r/ rc+1 n \ i , \ 1 n n -\ , / — rc+1 . n n+\\ t r\ n 77 . 

Pijkl( V eijk- V eijk) l h + l + l u V<xJ + ( a ePx t + PePx, ) l2 = Pijk F e + 

ffl = l 

m— 1 m— I •-' 



m=l m=\ 

+ A;Z{Z[(vi„+v; B )/2] 2 }- 

e=l m=l 
e=l * e=\ 

with parameters a e and fi e =\ — a e , e = 1,2,3 , by technology of 
global iterations 



n s,n+\ r / — s,n+\ , n s,n+\\ r\n i s,n+l 

P^eijk =Ha e P Xe +P e P Xe )-Q eij kK + V V eijk 



a=^'(21) 



e = 1,2,3,1 < j < AT, -1,1 < ; < N 2 -1,1 < it < N 3 -1, 
which assist in difference continuity equations in inward nodes 
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+[a 3 (p"vr i ) l3 +A( / )"vr i ) i3 ]=o, (22) 

i<* < iv x —1,1 < 7'<Ar 2 -l,l<fc<Ar 3 -l, 

Iterations in boundary nodes j = : 

+ (A- + io*0l« + io* - PmPuok ) 7 "im + (Ait v 2,-u - 

-/^t^V^+(/4 + i^Si + i-^'ot9 , £i) / ^t + i} = . (23) 
i =1,...,^, -U=l,.. .,iV 3 -l,j = 0, 

Iterations in boundary nodes & = : 

+ (A + l;0fli + l;0 " A/0 ftj/0 ) 7 ^i+1 + (Pff+10^+10 " 

"tfo ^V^+C^v^-^o <^)/^ 31 }=0, (24) 

i =l,...,N l -l,j = l,...,N 2 -l,k = 
Iterations in boundary nodes 1=0: 

, / n s,n+\ n n+1 \ / 7 , / n n+\ 

+ (Pijk v n jk -Po j k<PiOjk) /h ii+(Po j+ ik<P2 aj+ ik- (25) 

- PojkPlZly h 2j + l + (P0jk + l<P3W -P0jk<P£}jk) /h 3k + l} =0 * 

i = 0,j = l,...,N 2 -l,k=l,...,N i -\ 
Iterations in boundary nodes i = N l : 
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m+K^(p^r-p^)io+[p^ri(RT^-^ jk ]/T n+[+ 

+ \PNJk ( PlNJk FN^ljk V lN l -ljk)' 'Hn i + \FNJkY2NJk (26) 



.i+l \ / 1. i / -« ,.«+l ,.« n+1 

Vi* ~PN ] jk-mN ] j. 

i = N l ,j=l,...,N 2 -l,k = l,...,N 3 -l 



Iterations in boundary nodes j = N 2 : 

+ ( Atf 2 * PlW 2 * ~ A-ltf 2 t PlMtf 2 * ) ' K + (PiN 2 k (PliN.k ~ 

n 5, n+1 \ / / , / n n+i n n+i \ t i 1 a 

- PiN 2 -lk V 2iN 2 -lk ) ' h 2N 2 + (PiN 2 k<Pw 2 k ~ PiN 2 k-l<P3iN 2 k-l ) ' M = > 

j = N 2 ,l<k<N 3 -l,l<i<N l -l (27) 

Iterations in boundary nodes k = N 3 : 

m+^(pT-p^ ,0+ ^T l/ (^-p^^ + 

+ ( A//V 3 0W, ~ A-I./N3 ^1<-1.,7V 3 )/«li+ (PyiV 3 ^V 3 - (28) 



- PHJV, ^HJV, ) / ^2; + (Pjv, ^"p, - Af 3 -1 V 3? 3 -1 ) ' ^3 ^ = 0> 

l<i<iV 1 -l,fe = iV3,l<;<iV 2 -l 

By values « e = 0.5,/? e = 0.5,e = 1,2,3 from (21) - (28) prove 

scheme with global iterations (13) - (20). By other values of 

parameters a e and fi e = 1 - a e , e = 1,2,3 take place analogous 

schemes, which differ each other interconsisted approximations of 
gradients of pressure and continuity equation. 

In all without include applied above iteration processes, exactness of 
calculation of field of pressuse limited by exactness of execution of 
continuity equation. 
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Notice. Equality of zero some parameter a e ,/5 e , e = 1,2,3 means, 

that in this nodes continuity equation isn't applied for execution of 
pressure. 



§5. Technology building of scheme for march equations of 
compressible perfect gas 

March method of calculation of flow is applied in those cases, 
when in equations 2 derivatives by one from spatial coordinates 
possible disregard. For example, believe, that axle "x," direct by 

flow, put equal zero d v e /dx l = 0,e = l,2,3,3 T/dx 1 =0 or 
disregard diffusion 

d dv d dT 
flows — (// — -) = 0, (A ) = 0. In orthogonal to flow section 

dx i dx l dx l dx l 

Xj._j consider well-known (from boundary conditions or already 

calculated by algorithm) values desired functions V ej _ l j k ,T iA j k ,p iA j k 

and other, need find them values in section xi, that is v ei , k , T t - k , p ijk 
and other. For calculation stationary flows difference schemes 161 
apply in type iteration algorithm, when decision V eijk ,T jjk , p ijk and 
other, prove as limits 

lim v eijk =v em , lim T t J k =T ilk , lim p£ k = Pijk , 

n — >oo n — >oo n — >oo 

vj 1 -v n eijk Ta X ~Vk 
by execution conditions lim — — = 0, lim — — = 0, 

n— >co 7- n— >oo 7" 

and the like, it is enough to put in indicated schemes flows equal zero 

— (ju^) = 0,e = 1,2,3, — U— ) = 0, 

ox l ox l ox l ox l 
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and convection terms by direction x i approximate difference 
backwards 



dx l 



K "^ V = v m (v.* - C 1>k ) / K = v," v£ , 



that is put in schemes /6/ or" =0,5 + 0,5«',gn(vL ) = 1, as 
v i"# > 0' hk e manner (v"dT n I dx x ) ijk = v" ■ T" and so on. 

Constituent gradient of pressure dp I dx x 'm section x l5 approximate 

difference derivative backwards pl + , that is in schemes /6/necessary 

put cTj = 1, J3 X = , for two other constituent dp I dx m , m = 2,3 apply 

principle interconsisted approximation with continuity equation. In 

. . . ,d(pv,). ... . . 

continuity equation ( — ). ft approximate in section x 1( in view 

dx 1 

direction of flow in positive direction of axle xi also difference 
backwards: 



(_ a^ ) * = (pi ^ = (/ ^v - Pi-ijk v u-ijk)K\ 



where p t _i * v i,_i * calculated or well-known magnitudes, therefore 
in them don't put upper index n. In lineage march methods upper 

n 

index n don't put: V e i-\jk ~ V e i-ljk an d so on - 

On every iteration layout t n+l necessary execution inward 

n+\ .... 
iterations for executions p t - k as limit iteration 

lim» i .f +1 = »«+ 1 

additional upper index s conform to number of s - iteration. For 

n 

beginning field of iteration s = is selected value p& : 
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p:; +1 =pl,0<i<N v 0<j<N 2 ,0<k<N 3 . 

§6. Method of determinants in explicit iterationless scheme of 
solving stationary march equations 

Applied is stationary system of march equations 

m= i ox m dx l dx { 

m =2dx m dx m dx l 3 m=2 ox m 

^ dv e dp ^ d , dv 

m= i a* w a*. m=2 a* m a* m 

a 1 3 



? X^' )S fr )] ' e = 2 ' 3 ' 

m=l "* m m=2 t» m f-* m e=l m=l «* m 

1 , 

- pdivv -(—ju-ju')(divv) , 

3 



P^ 1 + ^ ^ + X ^ = 0, p = RpT, M = ju(T), X = MT) 

dx l dx l e=2 dx e 

where substitution is preliminarily made 

dp I dx 1 = RTdp I dxi + RpdT I dx l , 
Shortly denoted are backward difference derivatives 

u p =p Xl =(p iJk -Pi-^VK, 

approximating dT I dx x , dvyldx x , dpldx x . 

For calculations U T , C/j , U ' the following linear system 
comprising three equations was made, composed of scheme for v, 
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A-l* i( V U-ljk U l + E Vljl K m )i-ljk ] + ^-1 A + R PiAjk U T = 

m=2 

3 3 

= Pi-i jk F li -i jk + I^A W -X[(///3-//')v m , n ),],, 1# , 

ra=2 ra=2 

from schemes for continuity equation 

3 

VuAjkU p+ Pi-ljlUl+Y^eOi-ljk + (PVex e )i-lj k ] /2 =Q , 

e=2 

and energy balance 

3 3 

Pi.l jk [Vy_ ljk U T + Yj V mi-ljk ( T x m )i-ljk ] = H^m T x m )x m l-ljk + 
m=2 m=2 

3 3 



e=l m=2 

e=1 

-(V//') {(v,_ 1# -v,.,,)/^, + £[0^ +v Me ) M , /2]} 2 , 

■J c=2 

i > 1,1 < j < A^ 2 - 1,1 < >t < A^ 3 - 1 

In the initial cross-section it is assumed 

V e-ljk = V eOjk =<Pe0jk> e = 1 > 2 > 3 - 
Obtained was linear system of type 

a n U l + a l2 U + a 13 U T = b x , 

a 21 U 1 +a 22 U p +a 23 U T =b 2 

a 3l U 1 +a 32 U +a 33 U T =b 3 
Is solved by method of determinants, then found are unknown 
quantities in section x u for all j,k: v Ujk ,p ijk , T ijk : 

V = 1 W +h u l U vP, j k =Pi-i jk +h U P > T ijk = T u-ijk +K U t> after that, 
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by explicit schemes calculated are p ijk = p^RT^ and velocity 
components v e ijk » e = ^3 ; 



A* tO^* ( V ey * " V ei-ljk ) / *li + Yj V mi-ljk ( V «,„ )i-lj* 1 + 

m=2 

+ (P Xe +P x M' 2 = P ijk F m + 

3 3 

+ Z [ (a« v < ).v,„ w - X [ C" /3 - //>„*„, ), e ],-_!* * = 2 3 

m=2 m=l 

In initial conditions, at input, it is necessary to set density 

V l = V input ' V 2 = V input ' V 3 = V input ' "* = * input ' P = P input • Under 

flow periodicity, for instance, by x m ,m = 2 or m=3 index 
i m changed within limits 1 < i m < N m . In boundary nodes 
i m =0,i m = N m ,m = 2,3 calculation of density engages: either 

continuity equation, or condition of symmetry, on solid impermeable 
surfaces extrapolation formulae or other physical conditions. If in 
energy balance equation pdivv=0, there arises a system comprising 
only 2 equations 

a n U l +a 12 U p =b x , 

a 2l U x +a 22 U p =b 2 

Main disadvantage of march methods of such type is that for 
calculation of pressure are not applied difference analogues of §3 of 
canonical equation 

— y(- J — ) = > { — f- + — ^-[(-ju-ju')divv] + 

dt 2 RT tT dxf dxf 3^ 

d 

H [d/vC/TV; v ) - div(jugradv i ) - pF, ] } 

dx i 
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Chapter 11. UNIVERSAL METHOD OF ACCELERATION OF 
CONVERGENCE OF ITERATION ALGORITHMS 

Let for find decision p equation 

Ap = f, (1) 

where/- given quantity, A - nondegenerate linear operator, for 
which is inverse A~ , apply some coveregence iteration algorithm: 

B(p s+1 -p s )/r=Ap s -f,s = 0,l,..y (2) 

As a rule, with the exception of p s more part of proved information 

about decision p with previous steps of iteration p v ,v = 0,l,...,s-l 
rests without further applying in mind that don't use return. 

Basic idea offer universal method concludes exactly in return using 
of results p s ~' + of previous iterations (see Jakupov K.B. Ill) 

n n 

P S+2 = YjiP^'Yi = const ' & =1. n < s + 1, (3) 

as contains with some kind of degree of exactness of valuable 
information about desired solution p equation (1). 

In (3) coefficients y i = const are found from condition of 
minimization square of gilbert's form of residual 

ii n2 ii n2 

||fl' +2 || =JAp s+2 -f\\ D , 

some self-conjugate positive certain operator-generated D, for the 
sake of simplicity suppose in further D=E, E - identical (unit) 
operator. 

Definition. Acceleration, made by lineage combination (3), named 
m-cyclical acceleration, where m - quantity of different from zero of 
parameters y t , which are subject to calculate with indicated method. 
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This definition supposes opportunity of preliminary selective 
equalization to zero some parameters y- t ,i ^0. 

Let p s - ^-iteration/?; R" = Ap s — f - residual ^-iteration; 

/, p, p s — essence component of gilbert's space H; (f,tp) — scalar 

product in H; form of component \\p\\ = -Jip, p) ; 

H is field of definition and field of values of operator A. 
Theorem 1 (basic). In m=n+l — cyclical acceleration 

n n 

p s+2 = X YiP s ~ M ' y = const ' Z y ■ = l > l - n - s + 1 ' (4) 

i=0 (=0 

parameters y i ,i=\,...,n calculated as decision of system n lineage 
equations 

y\R s ~ M -R s+1 ( + Y J y i (R s ~ i+1 -R s+l ,R si+l -R s+1 ) = 

j=1 (5) 

j*i 

= (R s+1 ,R si+1 -R s+1 ), i = l,...,n, 
are necessary conditions of existence of extremum of function 



\\R s+2 \\ =\\Ap s+2 -f\\ = 



A^y.p^-f 



by this take place inequality II R s \\<\\ R s II , that is residuals 
don 't increase. 

Proving. Take action by operator A on two parts (4) and subtract 

f: 

Ap+ 2 -f= A±y iP ^ -f = ± 7i Ap s - M - f± 7i = 

i=0 i=0 (=0 

= I>,-(^-'' +I ~f),Ap s+2 -f =f j y i (Ap s - M -f) 

i=0 i=0 

Enter in last expression of residual, find 
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jr* =r R s+l +Zr i R s - M = (i-Y J r i )R s+1 +T J r i R s ~ i+1 = 



1=1 i=\ 

= R S+1 + J]y i (R s ~ i+1 -R s+1 ) 

!=1 

Further, by definition of form, have 

— ' ' II II 

=1 1=1 



\\R s+1 f =\\R s+1 f +J^ r f k s - i+1 -R s+1 f +2^ 7i . (R S+1 ,R S - +1 -R s+l ) + 



+^ZTr i r j (R s ~ M -R s+1 ,R sj+l -R s+1 ) (6) 

From expression (6) follow necessary conditions of existence of 

II s+211 2 
extremum d\\R I dy t =0,i = l,2,...,n in type of system of 

lineage equations (5). Apparently, enough conditions of minimum 

II l|2 

^21^5+2 I Qy* >fj,z =1,2,. ../i also are executed. Execution of 

inequality II R s+ \\<\\ R s+ His proved by method of mathematical 
induction below. Show more importance consequences from basic 
theorem. 

Consequense 1. In di-cyclic (m=2,i=l) acceleration 

p ,+2 =r p' +1 +riP'>r +ri=i (?) 

parameters y , y x are executed by formulas 

(R k+l ,R k+1 -R k ) , (R k ,R k+l -R k ) 

p*+l n« P* +l P* 

_ 

For conclusion of formulas (8) enough to put in (5) (5) n=l,i=l. 

Acceleration (7) equivalent of method of minimal residuals (m.m.r.) 
Krasnoselcky-Krein. In this method 

p' +l = p s + r s B~ l R s ,R s+l =R S + r s AB~ l R s , (9) 
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parameter r s is found from condition of minimum /? s+ and 
equals 



r = 



(R s ,AB~ l R s ) 



Wab-'ri 



(9') 



With other side, in view (7) and (9) 

p' +2 = r p' +1 +r l p'=Mp'+f,B- l R')+r l p' = p'+rof.B- 1 R', do) 

where r s ^ r s and in view (8) 



R s+l -R s 



n r s =- 



(R\(R s+l -R s )r s )_ (R - ," 



) 



i? I+1 -/? s 



R s+l -R s 



(R^AB^R 1 ) 

lAB-'R't 



Indicate y r s = r s , find, that in (9) and (10) residuals are 

minimazed equally. This fact is evidence of using formula (7) in 
methods of minimal residuals give higher rate, which agree with 
speed one's method. For acceleration of convergence of methods of 
minimal residuals, multi cyclic equations and dicyclic equations of 
following types are efficiency. 



Theorem 2. In di-cyclic accelerations (m=2) 

P s+1 =y Q P s+l +y l P s - ,+ \y Q +y l =l,i = l,2,...,s + l 

parameters /q,/, are calculated by formulas 



(11) 



_ (R s+1 ,R s+l -R s ~ i+l ) 

~ II R s+1 — R s ~' +1 II 2 _ 

by this take place inequality II R s \\<\\ R s II . 

P r o v i n g. In view (1 1) 

R s+2 =R s+1 + r .(R s - i+1 -R s+1 ), 
therefore, according with (12) find 



(12) 
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\\r s+2 ( = \\r s+1 + r .(R s - i+l -R s+1 )f = 



i, s . +1 |,2^ (R s+i ,R s,+1 -R s+i Y x 

= P(i-r — i^; — it). (13) 

\\K — K \\K 

By inequality of Koshi-Bynyakovskyi 

(R S+1 ,R S ~ M -R s+1 ) 2 <\\R s+1 f\\R s - M -R s+1 f, 
in view of (13) follow inequality II R s+2 \\<\\ R s+1 II . 
Theorem 3. In three-cyclic accelerations (m=3) 

P s+1 =/oP s+1 +LP S ~ M +7jP s - J+ \i* j*0,r +/i +/J =1 > (I 4 ) 
i, j — 1,2,. . . ,s + 1 parameters ^ , ^ , y. are executed by formulas 

(RS+lR s +i _^ i) |^ +i | 2 -(r s+ \r;:)j(r:%,r::) +1 ) 

Yi= A ' 

(RS+i RS+l _ R ^ ^ip _ (/? , +1 > ^ i)(/?;: i +] , ^ ) 

'' = A ' 

R::i +l =R s+i -R s - k+ \k=ij, (is) 

S+2|l II !-» A'-l-l II 



P r o v i n g. In view (14) 



^ +2 =^ +1 -7^-^ + i' 



.S+2II 2 



\\r m \\ =\\r s+l \\ +r-\\R:-L\\ +rj\\R s s%\\ -2 ri (R s+l ,R s ;4 +1 )- 

-2 rj (R s+ \R::) +l )+2y irj (R:%,R::) +l ) m 

Put (15) in (16) and make appropriate transformation, arrive to 
equality 
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ll n v + 2ii 2 n„ t+ ill 2 „ A z B-2ADC + C l E^ 

\R s+2 \ =\R \ (1 2 ), (17) 

A|7?' v+1 | 
where 

a = (r s+ \r;%),b = \]p? J+1 (,E = \\r s s %( 

c=(jrM£j +l) , d=(r:%,r::) +1 ) 

For that in (17) executed inequality 

\\R s+2 \\l3\R s+1 \\, (18) 

enough prove, that 

A 2 B-2ADC + C 2 E>0 (19) 

so A 2 B > 0,C 2 E > 0, that possible two cases: 1) ADC < 0, 

then inequality (19) executed by apparent way; 

2) ADC > 0, in this case execution (19) isn't apparently. Proving 
that by ADC > take place (19) rests upon that by inequality 
Koshi-Bynyakovskyi \d\ < vfi * \E , therefore, 

A 2 B - 2ADC + C 2 E > A 2 B - 2sign(AC)^B ■ Je + C 2 E = 
= (AVfi±cV£) 2 >0,etc. 
Theorem 4. In four-cyclic accelerations (m=4) 

p s+2 = 7oP s+l + hP s ~ M + rjP s - 1+l + r t p s - ,+1 J*j*t*o, 
y + ri+y } + y t = 1 JJ,t = i,2,-;S+i (20) 

parameters y Q , y i , y, , y t are calculated by formulas 
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Yi 



Q 


B 


C 


P 


E 


F 


T 


F 


M 



7 



A 


Q 


C 


B 


p 


F 


C 


T 


M 



7, 



A 


B Q 


B 


E P 


C 


F T 



A J A 

)S+l Z> ? +1 



A 

s-k+l 



(21) 



ro =i-r, -Yj -YrK-L =R s+l -R s - K+ ,k=i,j,t, 

ii n2 

A-\\K, ;i1 || ,B-{K s _ M ,K s _ j+l ),L -i.K s _ M ,K s _ t+l ), 



E= R 



s-i+l\\ ' 

v + 1 II 2 



5+1 D-S+l 



S-J+l ' V \s- J+l ' \s-f +1 ' 

q=(r s+ \r:%),p=(r s+ \r::L).t=(r" +, .r: + _; +1 ). 



A = 



7+1 

ABC 
B E F 
C F M 



by this ll/? v+2 ll<ll/? s+1 II. 

Formulas (21) follow from (5). Use method of mathematical 
induction possible put proving, that in theorem 4 and in basic 

theorem take place inequality \\R i+ II<II7? S+ II, as for n-1 and n=2 
this proved in theorems 2 and 3. 

Really, formula (20) can be written in type, which resembles three- 
cyclic (14): 



s+2 s-+l , .v-i'+l , — s 

P =YoP +YP +YP 



(22) 



YP 



li_ p <-J+i + _ 



Yj+Y, 



r ' P s -' + \Y = Yj+Y t 



Yj+Y 



If y. = ~Yt ' tnat (^O) possible write in type 
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P S+2 = VoP S+1 +hP S - M +/jP'>P S = P S ~ J+1 - p S -' +l (23) 
Formulas (22), (23) formally agree with (14), and that parameters 
Yo ' Yi ' Y' Yj are calculated from conditions of extremum. 

Theorem 3 is right for three-cyclic acceleration, therefore 
II R i+2 \\<\\ R s+1 II , etc. Similar arguments method of mathematical 
induction lets find, that in basic theorem will take place necessary 
inequality II R s+ \\<\\ R i+ II. Notable property of acceleration (3) 

is that inequality II R s+ \\<\\ R s+ 1 1 executed for any nonvacuous 
lineage operator A, from which not require so properties, as diagonal 
prevalence, self-adjointness, positive definiteness and etc., need only 
convergence of iterations. 

Numerical experiments, taken with two-cyclic acceleration by 
i > 5 showed, that in some iteration methods acceleration makes 
them convergence from 1,5 till 9 times faster. In the capacity of 
illustration below showed tables of exchanges of maximum-form of 
residual on example of decision of Dirichlet's task for Puasson 

equation on network < n <N V < m <N 2 ,N l = 20,N 2 = 20 with 
two-cyclic acceleration of method of Krasnoselskyi-Krein (m-K-K). 

In theorem 2 used only two previous iterations p s+ and p s ~ l+ , 
for example, for i=3 p s+ =Y p s+ +f 3 /"' for i=5 

p s+2 =YoP s+1 +Y 5 P s -\etc. 

For example, applying of two-cyclic acceleration for i=5 made by 
algorithm: 

1 2 3 4 5 -6 7 8 9 10 11 12 -13 s-4 j-3 s-2 s-1 s s+1 -s+2 

P ,p ,p ,p ,p ,p ,p ,p ,p ,p ,p ,p ,p ,...,p ,p ,p ,p ,p ,p ,p ,,.., 

Overlined iterations p ,p ,p ,...,p s+ are calculated by formula 
(11), but non-overlined 

1 2 4 5 7 8 9 10 11 12 s-4 j-3 s-2 s-1 s s+1 

P ,p ,p ,p ,p P ,p ,P P ,p ,...„P ,P ,P ,P ,P ,p ,..., 

calculated by basic (that is speed up) algorithm. 
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Experiments showed, that advantage by time takes place for two- 
cyclic accelerations with number i > 5 . 

Necessary take into account, that residuals in boundary nodes by 
edge conditions of Dirichlet, wittingly equal zero, by edge conditions 
of type fon Neiman calculated by given algorithm. Effect of 
acceleration is reached by calculations with two exactness. 

Tables of decrease of residual forms \\R S 11= max I R s nm I 

by two-cyclic acceleration p s+2 = y p s+1 +/iP' s ~' +1 

with exactness of iteration \\R S \\<s, £" = 0,01 



Iteratio 
n 

s=... 


(m-K-K) 
without 
acceleratio 
n 

\\R S \\=..., 
/=194 


Acceleration 

i=5 

\\R S \=..., 
s'=86 


Acceleratio 
n 

i=7 

\\R S \=..., 
/ =97 


Acceleratio 

n 

i=8 

\\R S \\=..., 
s* =104 


s=0 


806,820 


806,820 


806,820 


806,820 


s=10 


70,398 


73,002 


74,087 


74,903 


s=20 


55,096 


41,048 


45,618 


44,292 


s=30 


34,500 


14,071 


20,214 


11,872 


s=40 


20,966 


3,299 


6,903 


4,039 


s=50 


12,696 


1,067 


1,836 


1,934 


s=60 


7,685 


0,274 


0,828 


0,542 


s=70 


4,651 


0,047 


0,182 


0,297 


s=80 


2,815 


0,016 


0,074 


0,079 


s=90 


1,704 


0,005 


0,019 


0,031 


s=100 


1,031 


,'=86;0,001 


,-=9v;0,007 


0,011 


s=110 


0,624 






s" = 104 ; 


s=120 


0,378 








s=130 


0,229 
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s=140 


0,138 








s=150 


0,084 








s=160 


0,051 








s=170 


0,031 








s=180 


0,019 








s=190 


0,011 








s* =194 


0,006 
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Chapter 12. ABOUT APPLYING OF SCHEME OF 
PREDIKTOR-CORRECTOR IN HYDRODINAMICS 

In algorithms of alternative directions (as type Douglas-Rackford, 
Pisman-Rackford, Marchuk-Yanenko and others) as opposed to 
explicit and implicit schemes, arise problem of raising of boundary 
conditions for intermediate grid functions. This problem is worsen 
for boundary conditions as type fon Neiman, and becomes insolvable 
in fields with curved boundary. With this point of view effective is 
applying of conditionally stable explicit schemes or absolutely stable 
implicit schemes. 

For system of equations of hydrodynamics of viscous fluid, in view 
their irregularity, idea of applying of implicit schemes is realized in 
combination with semi-explicit schemes, technology of building of 
them is stated in Chapters 9,10. 

Technology of scheme of predictor-corrector is demonstrated for 
starting-boundary task: 



P( 



dt 



m=\ 



^ e)| dp 
dx m dx„ 



pF e + 
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+ I^(^>-f [(^-/Ot^l. -1,2,3, 

m= idx m 5x m 5x e 3 m=1 5x m 

,ar ^ dr. ^ a ,,&r. 

sr m= i 5x ffl m=1 ax m ax m 

at e=1 3x e 

v , I r=o =** e . v e | s =#> e > g = U3, T\ t=0 =d T ,T\ s =(Pr 

Predictor is monotonous scheme (see Chapter 6,9,10): 

Pijk V eijk = V-Px. ~ Qeijk K + l > e = 1,2,3, 

,P - /V I" 7 „ + l V # + Zji ~ V ei m + ~ V «„, J - F e J 

m=\ *- *- 



m=l m=l ~^ m=l 

-i i n i . « i n i n 

I v„,,,, I +v_, a I v_, a I -v_ 



(2) 



AW^s* -**)' V+2/ o T ^ + 1 t O* = 

m=\ *- *- 

m=\ m=\ 

3 3 3 i 3 

i . , n X ' X ' /, ." \2 „n X ' , .« / ■ .« ■ .'« \/X ' , .« \2 

+ /V Zj Zj (V - m ) ~ Piik Zj V *. -(^Mijk-Mijk) (X, V *e ) ' 
e=l m=l £=1 *^ e=l 

A/* ~^ +^(p"v" +1 ) e , g =0, 
1 < i < N, -1,1 < 7 < Af 2 -1,1 < k < N 3 -1, 
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k+1 _ n 

Pijk Pijk 



"n+l 



n+1 _ n 

Pijk Pijk 



n+l 



Hp n v n+ X t +J»" +1 U, =0,i e =0,l<i q <N m -1, 



m-\ 



+ (pV% +5>W =0,i e =N p l<i q <N m -1, 



m=l 



e = 1,2,3, 4 = 1,2,3, e * q, p n+l = Rp n+l T" +l , 



n+l 



v.. 



s =<Pt , 



M = M X -Mi = M y -Mi =M Z -A =A x >*2 =A y -A = A - 

Realization of so schemes is made by global iterations, that let in 
Chapter 9 and 10. In result calculated fields of components of 

n+l i a o n+l 

velocity vector and pressure V e ,e — l,Z,j, p on time moment 

K+i ■ 

Corrector is the same monotonous, but implicit scheme on time 

layout t n+2 . Corrector can consist from two types. 
Lineage corrector. It takes type: 

^ji+l n+2 ^Ji+l n+l r n+l . /~<n+2-i . _ „ ,_, 

Pijk v ei jk =P ijk v eijk ~[p Xe +G eijk ]T n+2 , e = 1,2,3, (3) 



Gk+2 n+l _R~'r 

eijk =Pijk Y^l 



3 I n+l i . n+l i n+l i n+\ 

i ' V ■; \+V ••, V •-, —V , -^ 

mijk mijk n+l mijk mijk n+2 i 77 I 

V- H V — r -1 

~ ex m r\ ex m J e J 



m=l 



n+l 



E/ n+l n+2\ X -1 / n+l \ n+2 , r/M \n+l \\^ n+l -1 

(M V ex m )x m -L(Mn, -M m JV eVm +[(— // )2j„u m X> 

m=\ m=[ ~* m=\ 

1 < i < N, -1,1 < ; < N 2 -1,1 < k < N 3 -1, 



n+2 



v„ 



n+2 1 /-> o 

=#, ,e = 1,2,3 
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In this scheme well-known grid magnitudes are p n , p , T" , 

ju" = ju(T" ) and components of speed V e , e = 1,2,3 , they are 

found by implicit scheme (2). In (3) difference derivatives taken on 

moment of time t n+2 , but coefficients by them on layout are t n+l . 

In this case equations (3) are system of lineage algebraic equations, 
concerning components of speed 

v n e ^,i=i,...,N l -ij=i,...,N 2 -i,k=i,...,N 3 -i 

on time's layout t n+2 . Though system (3) has diagonal prevalence 

«+l n+2 

because of presence term of type ^— , that is executes 

T n+2 

condition of convergence of iteration of Jakobi method. 
Present (3) in compact type 

+ P?+G£=0,e = \,2& (4 ) 



n+i / n+z n+i \ 

Pijk \ V eijk ~ V eijk > t ^n+l , /-.n+2 



eijk 
v n+2 



l<i<N l -l,l<j<N 2 -l,l<k<N 3 - 1, 



n+2 



S, 



=(p" ,e = 1,2,3 



Well, s- iteration indicated V ej - k " , in limit for convergence iteration 

I* /i + 2 S W+l 

process of Jakobi lUTt V eijk ' = V eijk . Starting field of iteration by 

n -. n+2,Q n+1 100. 

s = equals V e =V e ,e = 1,2,3. \ 
Calculate residuals s - iteration R s e , e = 1,2,3 : 



n+l/n+z,s _,.n+l\ 
D s Pijk V eijk eijk ' n+1 . /~<n+2,s 

K eijk = + P~x e + U eijk ' 



-n+2 



n+2, s 
1 A ' A — J — L " 2 A ' A — "-— iT 3 ^ y eijk 



\<i<N l -1,1 <j<N 2 -1,1 < ifc < N 3 -1, v, 



n+2 
- (I) 



S„ ~ Veijk 
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R 



eijk 



= V n+2 > S 
S h eijk 



n+2 
S h Veijk 



-<p:£=0,e = l,2,3,s=0X...,s* 



Extract coefficients, state by V ei j k in G" +2 

/ vw \ , / n+l n+l \ I r\ i n+l i r 

„+i (/*, -M m J + (Mi + i jk +Mi jk ) /2 , '«» !-^ 



n+l I n+l 



-+!// 



^n+1 



fcrfAi 



- + /?- 



p xi+\ 
/ vw \ , / n+l , n+l \ / o I n+l i , n+l 

(M x -M m J + (Miik +A-i*) /2 l«« \+ u a 



Kfi* 



+ p- 



2*, 

r+1 I 



2/i, 



L +2,u ~~ 



/ vvm \ , / n+l , n+l\ i r\ i n+l i 

(^ -MmJ + iUj+U+Mijk ) /2 , |V ^ '-'SI 



n+l | ..n+l 



^A 



+ /> 



2/j.. 



n+l I , n+l 



/ wu \ . / n+l . n+l \ / ^ i n+l i , 

c „+i = (M y -M m J + (Mij k +Mij-ik) /2 { \ v uk \+ v i jk 



-2// 



K h yj 



2h 



yj 

+1 I . . n+l 



/My \ , / n+l , n+l \ i r\ , n+l i 



C +3,u — 



W** 



+ P- 



afc+1 



n+l _ v >~z /~min 



M* 



+p- 



2^ 

Z +1 1 



2ft 



'afc 



,«+2,s+l 



In Jakobi method s + 1 - iteration V c " ,e = 1,2,3 is identified by 
algorithm: 



o n+l 3 

,.n+2,.s+l ,.n+2,s ns /r^!/' 1 , X -1 ( „n+\ , „n+l \1 

^"n+2 m=l 

f = l,...,iV 1 -l,j=l,...,iV 2 -l,fe=l,...,iV3-l, 



n+2,i+l 
eijk 



n+2 



= (Peijk > e = 1» 2 »3. * = 0,1, . . . ,5 
Iterations are stopped by execution of criterion 



max I R s . 1k l< s, s « 0, £• ^ 0, e = 1,2,3 



l<y<7V 2 -l 



<?y& 
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Interchange is made by this way: on layout of time t n+l realized 

predictor, on layout t n+2 corrector. Of course, possible apply other 

iteration methods. For example, comfortably to write method 
Libman-Zeidel with help next difference operations: 

A n : 1 =c; + jA +m -c n _ + m 1 A_ m ,m = 1,2,3, where A +m =T +m -E,m = 1,2,3 - 
operator of formation of differences forward, A_ m =E-T_ m ,m = 1,2,3 
- operator of formation of differences backward, E - unit operator, 
T +m ,m = 1,2,3- operators of shear forward, T_ m ,m = 1,2,3 - 
operators of shear backward. 

System (4) is written in operator type 



■crvv" -v" +1 ) 3 

' v e e \ 1 / n+l a n+1 a \ n+2 n+l . n+l r? 1 1 o 
= Z>+mA™" C -™/->* ~i\ +P ^>= 1,2,3, 



'n+2 m=l 



for which method of Libman-Zeidel is realized as 

n+l/ n+2,s+l w 



_n+l/..n+2,j+l ,.n+l\ 3 



r n+2 m=l 



Y\ n+l /t n+2,.s n+2,s+l\ n+l 1 n+2,.s+h n+l , n+l 77 



i = l,...,N 1 -l,j = l,...,N 2 -l,k = l,...,N 3 -l, 



n+2,s+l 
e 



S h 



= ^" + V = l,2,3,s = 0,l,...,s 



Similar way can use method of upper relaxation and other methods. 
Particularly, considerable effect gives applying of universal method 
of acceleration of convergence iterations Chapter 11. 

Temperature T n+ on moment of time t n+2 also calculate by 

method of Jakobi (or other method) as decision of next system of 
lineage algebraic equations with diagonal prevalence: 

3 I v "+ 2 I +v "+ 2 I v "+ 2 I _ v "+ 2 
„"+ 1 „ r/T" I + 2 Tin+U I _ , V" 1 / m 'i k m 'J k t"+ 2 i m 'J k m 'J k t"+ 2 m 
Pijk C v\-( T ijk - T ijk V^+z+Z/ o ^ + ~ ^ )] = 
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=£a"'C\+Ear 1 -4,„>C: + 

m=\ m=\ 

3 3 3 1 3 

i .."+ 1 X~''W,,n+2\2 „"+lV„"+ 2 ( . ."+ 1 ,,'»+h/'V,.' I + 2 \ 2 

+ /V 2j2j (V «J ~^fl* LP&. -(-Mijk -Mijk )(2 J V - e ^ 
e=\ m=\ e=\ ^ e=l 

Calculate residual s - iteration R^ : 

^ = £fe C v [(^ - ^ ) / T n+2 + 

3 I v "+ 2 I +v "+ 2 I v "+ 2 I _ v "+ 2 

. V^ / mi/t mp jm+2,s , mijk mijk t»+2,s\-i_ 

m=\ ^ ^ 

-t^xTx + ix 2 -^ki: -<iix 2 ) 2 + 

m=\ m=\ e=\ m=\ 

In method of Jakobi s + 1 - iteration T " is determined 
by algorithm: 

«+l 3 

7i«+2,s+l rpn+2,5 r>s irrijk v V -1 /n+2 , „n+2 \-i 
^ = T ijk ~ R eijk /[ + 2>W + C ~^ )1 ' 

7 n+2 m=l 

i =1,...,^, -l,j = l,...,N 2 -l,k=l,...,N 3 -1, 

,_n+2 ai * 

s ft ~~ 'Prijk » * ~~ u » i » • • • > s 



rjin+2,s+\ 
I ijk 



n+2 

hjk ' 

n+2 n+2 



There are coefficients C+ m x » ^-ml are similar by coefficients 

n+2 n+2 ^"+2 ~w+2 

C + , C_ . In lineage corrector in c + m x ■> c - m A use d coefficient of 
thermal conduction A" +1 =A(T n+l ) on layout t n+1 . 
Iterations are stopped by execution of criterion 
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max I R^jk \^ £, £ ~ 0, s ^ 

l<j<N 2 -l 

l<k<N 3 -l 

Non-lineage corrector. It has type: 

n+1 n+2 _n+l n+1 r 1+1 ■ /~^«+2-i_ , ,. _ 

P V e =P V e -\-P~r + G e K+2> ^ = 1,2,3, 



g: +2 =p^[ 



3 , v „ +2|+v „ + 2 



»+2 IV„ I V m v » + 2-,_ F ^ 



,n+l 



-5> n X\ -lwr 2 -/u.)<i. +[(^-- ^ +1 )I^;l v 



/»-i 



m=l *^ m=l 

1 < i < N, -1,1 < ; < N 2 -1,1 < k < N 3 -1, 



v.. 



n+2 



n+2 1^0 

=#, ,e =1,2,3 



,n+2 



In // m ,m = 1,2,3 used coefficient of viscosity fj" 
components of speed V e , e = 1,2,3 . 

Execute residuals s - iteration R e , e = 1,2,3 : 



and 



^eijk ~ 



n+1 / n+2,i n+l\ 

M/& v'gjji ep / , „»+l , p^n+2,,s 



+pr+G, 



eijk 



''n+2 



Ui^-l.U./'^-U^Jfc <#,-!, 



n+2,i 



n+2 
S h ~^eijk • 



R* 

eijk 



= V 



n+2,s 
eijk 



<p:£=0,e = l,2,3,s=0,l,...,s* 



n+2,s+\ 



In non-lineage corrector s + 1 - iteration V e " , e = 1,2,3 is 
determined by similar algorithm of Jakobi: 
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«+l 3 

= V eijk ~ R eijk /[ + lS C +W + C -M )] ' 



n+2,j+l _ ,,n+2,.s D .5 /i^ijk 



ei J k ■ • _ 

t n+2 m=l 

f = l,...,iV 1 -l,j = l,...,iV 2 -l,it=l,...,iV3-l, 



n+2,.v+l 
eijk 



Sh = (p"y k 2 ,e = l,2,3,s =0,l,...,s* 
Iterations are stopped by execution of criterion 

max I R s * ijk \<s, s^O,s^O,e = 1,2,3 

l<j<N 2 -l 
l<k<N 3 -l 

Algorithm of calculation of temperature T n+ on moment of time 
t n+2 is realized after definition V e ,e = 1,2,3 and indicated above 
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The upper figure shows vortex separation inflow around buildings 
with different height. Uniform stream runs from the left side. Grid 140x80. 

The lower drawing presents a picture of plate cross flow. 
Developed behind the obstacle two vortexes develop with time into 
interacting with each other and flow -drifted nonsymmetric vortexes. 
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Re=40 000, 140x80, t=9685*tau, tau=0.001 



>" 1 




Field of actual velocity vector V = V + V in turbulent flow around four 
buildings, top view. Dimensionless averaging time is equal Dg=2*T, 
T =tau. 
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